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Algorithms for the computation of the forward and inverse geodesic problems for an ellipsoid of revolution are 
derived. These are accurate to better than 15 nm when applied to the terrestrial ellipsoids. The solutions of 
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1. INTRODUCTION 

A geodesic is the natural "straight Une", defined as the 
l ine of mi nimum curvature, for the surface of the earth 
jHilbert and Cohn-Vossenl Il952l pp. 220-222). Geodesies 
are also of interest because the shortest path between two 
points on the earth is always a geodesic (although the con- 
verse is not necessarily true). In most terrestrial applications, 
the earth is taking to be an ellipsoid of revolution and I adopt 
this model in this paper 

Consider two points A, at latitude and longitude {(j>i, Ai), 
and B, at {(t)2,^2), on the surface of the earth connected by 




FIG. 1 The ellipsoidal triangle NAB. N is the north pole, NA and 
NB are meridians, and AB is a geodesic of length S12. The longi- 
tude of B relative to A is A12; the latitudes of A and B are (l>i and 
(j)2- EFH is the equator with E also lying on the extension of the 
geodesic AB; and ao, qi, and Q2 are the azimuths of the geodesic 
at E, A, and B. 
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a geodesic. Denote the bearings of the geodesic (measured 
clockwise from north) at A and Bhy ai and a2, respectively, 
and the length of the geodesic by S12; see Fig.IU There are two 
main "geodesic problems": the li/recf problem, given 0i, S12, 
and ai determine (f)2, A12 = A2 — Ai, and a2; and the inverse 
problem, given the 0i, 02, and A12, determine S12, ai, and 012. 
Considering the ellipsoidal triangle, NAB, in Fig.[Tl where N 
is the north pole, it is clear that both problems are equivalent 
to solving the triangle given two sides and the included angle. 
In the first half of this paper (Secs.|2]-|9), I present solutions to 
these problems: the accuracy of the solutions is limited only 
by the precision of the number system of the computer; the 
direct solution is non-iterative; the inverse solution is iterative 
but always converges in a few iterations. In the second half 
(Secs. [T0l - [T5] l, I discuss several applications of geodesies. 

This paper has had a rather elephantine gestation. My ini- 
tial work in this area grew out of a dissatisfaction with the 
wi dely used a lgorith ms for the main geodesic problems given 
bv IVincentvl (Il975ah . These have two flaws: firstly, the al- 
gorithms are given to a fixed order in the flattening of the 
ellipsoid thereby limiting their accuracy; more seriously, the 
algorithm for the inverse problem fails in the case of nearly 
antipoda l points. Startin g with the overview of the problem 
given bv lWimamsl(l2002h . I cured the defects noted above and 
included the algorithms in GeographicLib (' Karnevil2010l) in 
March 2009. At the same time, I started writing this paper and 
in the course of t his I came across references in, for example, 
IRainsfordI d 1 955b to work by Euler, Legendre, and Bessel. Be- 
cause no specific citations were given, I initially ignored these 
references. Howeve r, when I final ly stumbled across Bessel's 
paper on geodesies (lBesselLlT825l) . I was "like some watcher 
of the skies when a new planet swims into his ken". Overlook- 
ing minor quirks of notation and the use of logarithms for nu- 
merical calculations, Bessel gives a formulation and solution 
of the direct geodesic problem which is as clear, as concise, 
and as modern as any I have read. However, I was surprised 
to discover that Bessel derived series expansions for the geo- 
desic integrals which are more economical than those used 
in the English-language literature of the 20th century. This 
prompted me to undertake a systematic search for the other 
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original papers on geod esies, the fruits of which are available 
on-line karnevl l2009h . These contained other little known 
results — ^probably the most important of which is the concept 
of the reduced length (Sect. |3]l — which have been incorpo- 
rated into the geodesic classes in GeographicLib. 

Some authors define "spheroid" as an ellipsoid of revolu- 
tion. In this paper, 1 use the term in its more general sense, 
as an approximately spherical figure. Although this paper is 
principally concerned the earth modeled as an ellipsoid of rev- 
olution, there are two sections where the analysis is more gen- 
eral: (1) in the development of the auxiliary sphere. Sect. |2] 
and Appendix lAl which applies to a spheroid of revolution; 
(2) in the generalization of the gnomonic projection, Sect.fTSl 
which applies to a general spheroid. 




2. AUXILIARY SPHERE 

The study of geodesies on an ellipsoid of revolution was 
pursued by many authors in th e 18th and 19th centuries. Th e 
important early pa pers a re bv IClairautI dlTSSl) . lEulerl dl755b 
[Dionis du Sejoi^ d 17891 Book 1, Chaps. 1-3), Legendrd 
(fl789t 11806)'. and'Orianij (1 18061 1 18081 II8IOI) . IClairautl(ll735h 
found an invariant for a geodesic (a consequence of the ro- 
tational symmetry of the ellipsoid); this reduces the equa- 
tions for the geodesic to quadrature. Subsequently. 'Legen drd 
([I8O6) and Oriani (1806) reduced the spheroidal triangle in 
Fig. [1] i nto an equivalent triangle on the "auxiliary" sphere. 
iBessell d 1 825^ provided a method (using tables that he sup- 
pUed) to compute the necessary integrals and allowed the di- 
rect problem to be solved with an accuracy of a few centime- 
ters. In this section, 1 summarize this formulation of geo- 
desies; more details are given in Appendix lA] in which the 
derivation of the auxiliary sphere is given. 

1 consider an ellipsoid of revolution with equatorial radius 
a, and polar semi-axis b, flattening /, third flattening n, ec- 
centricity e, and second eccentricity e' given by 



/ 

n 

„'2 



(a-6)/a, (1) 

(a-6)/(a + 6) = //(2-/), (2) 

ia^~b^)/a'^f{2-f), (3) 

(a2 - b^)/b^ = eV(l - e^). (4) 



In this paper, I am primarily concerned with oblate ellipsoids 
(a > 6); however, with a few exceptions, the formulas apply 
to prolate ellipsoids merely by allowing / < and < 0. 
(Appendix |D] addresses the modifications necessary to treat 
prolate ellipsoids in more detail.) Most of the examples in this 
paper use the WGS84 ellipsoid for which a = 6 378.137 km 
and / = 1/298.257 223 563. (In the illustrative examples, 
numbers given in boldface are exact. The other numbers are 
obtained by rounding the exact result to the given number of 
places.) The surface of the ellipsoid is characterized by its 
meridional and transverse radii of curvature. 



(5) 
(6) 
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FIG. 2 The elementary ellipsoidal triangle NEP mapped to the aux- 
iliary sphere. NE and NPG are meridians; EG is the equator; 
and EP is the geodesic. The corresponding ellipsoidal variables are 
shown in parentheses. 



respectively, where 



w 



(7) 



Consider a geodesic which intersects the equator, = 0, in 
the northwards direction with azimuth (measured clockwise 
from north) ao € [— ^tt, ^tt]. I denote this equatorial cross- 
ing point E and this is taken as the origin for the longitude A 
and for measuring (signed) displacements s along the geode- 
sic. Because this definition of longitude depends on the geo- 
desic, longitude differences must be computed for points on 
the same geodesic. Consider now a point P with latitude (j>, 
longitude A, a displacement s along the geodesic and form the 
ellipsoidal triangle NEP where N represents the north pole. 
The (forward) azimuth of the geodesic at P is a (i.e., the angle 
NPE is TT - a). 

Appendix lAl shows how this triangle may be transferred to 
the auxiliary sphere where the latitude on the sphere is the 
reduced latitude /3, given by 



tan (3 ~ (1 — /) tan (p 



(8) 



( iLegendrd,ll806l p. 136), azimuths (ao and a) are conserved, 
the longitude is denoted by uj, and EP is a portion of a 
great circle with arc length a; see Fig. m (ICavlevI (1 18701 
p. 331) suggested the term parametric latitude for (3 because 
this is the angle most commonly used when the meridian el- 
lipse is written in parametric form. ) Applying Napier's rules 
of circular parts (iTodhunterl Il87l[ §66) to the right triangle 
EPG in Fig.|2]gives a set of relations that apply cyclically to 



[ao, ^TT - cr, ^TT - a,/3,a;]. 



sin ao = sin a cos j3 
= tan UJ cot cr, 
cos a = cos (3 cos uj 
= tanao cot a, 



(9) 
(10) 

(11) 
(12) 
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cos a = cos w cos ao 
= cot (T tan /3, 

sin 13 — cos ao sin a 
= cot a tanw, 

sin uj = sin a sin a 
= tan/Jtanag. 



(13) 
(14) 
(15) 
(16) 
(17) 
(18) 



In solving the main geodesic problems, 1 let P stand for A or 
B with the quantities (3, a, a, uj, s, and A acquiring a subscript 
1 or 2. 1 also define ai2 — <^2 ^ <^i, the increase of a along 
AB, with ijji2, si2, and A12 defined similarly. Equation (|9]l is 
Clairaut's characteristic equation for the geodesic, Eq. (IA2]l . 
The e llipsoidal quantities s and A are then given by (Be sseU, 
[m Eqs. (5)) 



1 ds _ dA _ 
a dcr AljJ 

where w, Eq. O, is now given in terms of (3 by 



■ cos^ /3, 



(19) 



(20) 



and where the derivatives are taken holding ao fixed (see Ap- 
pendix|A]l. Integrating the equation for s and substituting for 
(3 from Eq. (115b gives (lBesselLll825l §5) 



^ = / Vl + Zc^ sin2 cr' dcr', 



where 



k = e' cos ao- 



(21) 



(22) 



As iLegendrj dlSl ll §127) points out, the expression for s 
given by Eq. (1211 1 is the same as that for the length along the 
perimeter of an ellipse with semi-axes b and by/1 + fc^. The 
equation for A may also be expressed as an integral in cr by 
using the second of Eqs. (IA4b . du/da = sin a/ cos/3; this 
gives (lBesse]Lll825L §9) 



X = io — f sin ao 



/ 



i + (i-/)vr+^ 



' sin^ <7' 



:da'. 



(23) 

Consider a geodesic on the auxiliary sphere completely encir- 
cling the sphere. On the ellipsoid, the end point B satisfies 
(j)2 — and a2 — cei; however, from Eq. ( |23T l, the longitude 
difference A12 falls short of 2tt by approximately 27r/sinao. 
As a consequence, geodesies on ellipsoids (as distinct from 
spheres) are not, in general, closed. 

In principle, the auxiliary sphere and Eqs. (|211 and (|2Jt 
enable the solution of all geodesic problems on an ellipsoid. 
However, the efficient solution of the inverse problem requires 
knowledge of how neighboring geodesies behave. This is ex- 
amined in the next section. 



3. REDUCED LENGTH 



Following Bessel's paper, iGaussI (1 19021) studied the prop- 
erties of geodesies on general surfaces. Consider all the geo- 
desies emanating from a point A. Define a geodesic circle 




FIG. 3 Geodesies from a point <^i — —30°. The east-going geo- 
desies with azimuths ai which are multiples of 10° are shown as 
heavy lines. The spherical arc length of the geodesies is (J12 = 180°. 
The geodesies are viewed from a distant point over the equator at 
A — Ai = 90°. The light lines show equally spaced geodesic circles. 
The flattening of the ellipsoid was taken to be / = i for the purposes 
of this figure. 



centered at A to be the locus of points a fixed (geodesic) dis- 
tance S12 from A; this is a straightforwa rd ext ension of the 
definition of a circle on a plane. iGaussI d 19021 §15) proved 
that geodesic circles intersect the geodesies at right angles; 
see Fig. [3] 

(Circles on a plane also have a second property: they are the 
curves which enclose the maximum area for a given perimeter 
On an ellipsoid s u ch curves h ave co nstant geodesic curvature 
(iMindind, Il830l) . iDarbouxl (1 1894 §652) adopts this as his 
definition of the geodesic circle; however, in general, these 
curves are different from the geodesic circles as defined in the 
previous paragraph.) 

feauss (1902) also introduced the concept of the reduced 
length mi2 for the geod esic which i s also the subject of a 
detailed investigation by IChristoffell (Il910l) (the term is his 
coinage). Consider two geodesies of length S12 departing 
from A at azimuths ai and ai +dai. On a flat surface the end 
points are separated by S12 dai (in the limit dai — > 0). On 
a curved surface the separatio n is mi2 dai where mi2 is the 
reduced length, feaussi (1 19021 §19) showed that the reduced 
length satisfies the differential equation 



dm 
d^ 



K{s)m = 0, 



(24) 



where K{s) is the Gaussian curvature of the surface. Let 
m{s; si) be the solution to Eq. (l24l i subject to the initial con- 
ditions 



m{si;si) = 0, 



d?7i(s; si) 



ds 



1; 



then mi2 is given by TO12 = 'rn{s2; si). Equation ( |24] | obeys 
a simple reciprocity relation m(si; S2) = —m{s2; si) which 
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gives the result ( IChristoffell[l910i §9) that the reduced length 
is invariant under interchange of the end points, i.e., 77121 = 
-mi2. 

For a geodesic on an ellipsoid of revolution, the Gaussian 
curvature is given by 

1 _ &2 ^ 1 ^ 

pv a'^w'^ 52(1 + A:^ sin^ cr)2 ■ 



iHelmertldlSSOl §6.5) solves Eq. (EDi in this case to give 
^12/6= Vl + /c^ sin^ (72 cos (Ti sin (72 



— Vl + fc^ sin^ (71 sin cri cos (72 

— COSCTi C0SCT2(J((72) — J((7i)), 



(26) 



where 



J((7) 



__fc2sinV^ 

Vl + A:^ sin^ cr' 
1 



:d(7' 



+ fc^ si 



:d(7'. 



(27) 



sm a 



In the spherical limit, Eq. ( 126b reduces to 

mi2 = asin(7i2 = asin(si2/a). 

iGaussI (Il902i §23) also introduced what I call the geode- 
sic scale M12, which gives the separation of close, initially 
parallel, geodesies, and which is given by another solution of 
Eq. (I24I 1. M{s\ si), with the initial conditions 

AM{s;si) 



M(si;si) = 1, 



ds 



= 0. 



iDarbouxl (Il894 §633) shows how to construct the solution 
given two independent solutions of Eq. (l24b which are given, 
for example, by Eq. (|26] | with two different starting points. 
This gives 



M12 



cos (7i cos (72 



sm (72 



: sm (7i sm CT2 



a/1 + fc2 sin^ Ti 

SinCTl cos (72 ( J((72) — J((7l)) 



(28) 



sm (7i 



where i\fi2 = A/(s2;si). Note that M12 is not symmetric 
under interchange of the end points. In the spherical limit, 
Eq. ( |28] ) reduces to 

M12 = C0S(7i2 = C0s(si2/(l). 

By direct differentiation, it is easy to show that the Wronskian 
for the two solutions mi2 and M12 is a constant. Substituting 
the initial conditions then gives 



M12- 



dm 



12 



dS2 



TO12- 



dMi2 
ds2 



1. 



(29) 



There is little mention of the reduced length in the geodetic 
literature in the English language of the last century. An ex- 



ception is iTobevI d 19281 Prop. IV) who d erives an exp ression 



. pr 

for TO 12 as a series valid for small 512/0 (lRapij.ll99lL §4.22) 



4. PROPERTIES OF THE INTEGRALS 

The solution of the main geodesic problems requires the 
evaluation of the three integrals appearing in Eqs. (|211 . (|27T i, 
and (|23] |. In order to approach the evaluation systematically, I 
write these integrals as 



h{cf) = / Vl + fc2 sin^ (7' da', 

JO 



Vl + fc2sin (7 



2 (7' 

2-/ 



d(7', 



sin^ (7' 



:d(7'. 



In terms of Ij{cr), Eqs. fTH . dZTl i. and ( l23T l become 

s/b^h{a), 

J((7) =/l((7)-/2((7), 

X = uj- /sinao /sIct)- 



(30) 
(31) 
(32) 



(33) 
(34) 
(35) 



The integr als Ijicr) may be expressed in terms of e lliptic 
functions as dForsvthl \ma IJacobiL flSSl iLutheil fl856l) 



(■Ml 

/i((7) = fcW nd^(w',fci)du', 
Jo 



/ sin^ ao Jo 1 + cot^ ao cd^ (m' , fci ) 
tan"""^ (sin ao tan a) 
f sin ao 



(36) 
(37) 

(38) 



where fci = fc/^/TTI^, fc^ = Vl - Af = 1/VTTF, 

am(ui — if (fci), fci) = (7 — iyr, 

cdfrr, fc) and nd (x, k) are Jacobian elliptic functions 
dOlverefaZl 2010l jS22.2) am(a;, A:) is Jacobi's ampHtude 
function dOlver et a/.Ll2O10l §22.16(i) ), and K(k) is the com- 
plete elliptic integral of the first kind (lOlver ef a/.ll2010l §19j 
f 2(ii)). The integ rals can also be written in closed form as 
Ce^ndri [Unr § § 1 27- 1 28) 



ci, 



Hi 

h{<j)^k[F{<J~^TT,kl)~C2 



where 

G{(j), a^,k) 



A 2 ^ \^ 2 
I sm ao 

tan^ ^ (sin ao tan a) 

f sin ao 

Vl - fc2 sin^ e 



G{a — — cot^ ao, fci) 



C3, 



(39) 
(40) 

(41) 



= 1- 



1 9-2 

1 — sm 

fc2 



d6' 



n(,/),a2,fc) + ^i^^((/),fc), (42) 
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the integration constants cj are given by the condition Ij (0) = 
0, and F{(f),k), E{(t),k), and 11(0, are Legendre's in- 
complete ell i ptic in tegrals of the first, second, and third kinds 
(lOlveref adBOlol §19.2(ii)). 

There are several ways that the integrals may be com- 
puted. One poss i bility is merely to utilize standard algorithms 
(iBulirschl 1 19651: ICarlsonl 1 1995b for elUptic functions and in- 
tegrals for the evaluation of Eqs. (|39| -(|4l ) . Alternatively, 
some authors, for example Saito ( 1970ill979l) . have employed 
numerical quadrature on Eqs. (|3Ut and (|32t . However, the 
presence of small parameters in the integrals also allow the 
integ rals to be exp ressed in terms of rapidly converging se- 
ries. Bessell (1 1825b used this approach and tabulated the co- 
efficients appearing in these series thereby allowing the direct 
geodesic problem to be solved easily and with an accuracy of 
about 8 decimal digits. I also use this technique because it 
allows the integrals to be evaluated efficiently and accurately. 

Before carrying out the series expansions, it is useful to es- 
tablished general properties of the integrals. As functions of 
a, the integrands are all even, periodic with period tt, posi- 
tive, and of the form 1 + 0{f). (Note that = 0{f).) The 
integrals Ij [a) can therefore be expressed as 

/, (a) = A, {a + (a)) , for j - 1, 2, 3, (43) 

where the constant Aj = 1 + 0{f) and Bj{a) = 0(f) is odd 
and periodic with period tt and so may be written as 



1=1 



Cjism2l,a forj = 1,2,3. 



(44) 



In addition, it is easy to show that Cji ~ 0{f^). In order to 
obtain results for s. A, and accurate to order truncate 
the sum in Eq. (|44] | at I = L for j = 1 and 2 and at I ~ 
L — 1 for j = 3. (In the equation for A, Eq. ( |35] |. I3{(j) is 
multiplied by /; so it is only necessary to compute this integral 
to order /^~^.) Similarly, the expansions for Aj and Cji may 
be truncated at order for j = 1 and 2 and at order /^~^ 
for j = 3. 

The form of the trigonometric expansion, Eqs. (l4Jt and 
(|44] |. and the subsequent expansion of the coefficients as Tay- 
lor series in / (or an equivalent small parameter), that I de- 
tail in the next section, provide expansions for the integrals 
which, with a modest number of terms, are valid for arbitrar- 
ily long geodesies. This is to be distinguishe d from a num ber 
of approximate methods for short geodesies (lRapp[[l99l[ §6), 
which were derived as an aid to computing by hand. 



5. SERIES EXPANSIONS OF THE INTEGRALS 

Finding explicit expressions for Aj and Cji is simply matter 
of expanding the integrands for small k and / ena bhng the 
integrals to be evaluated. I used the algebra system i Maximal 
(f2009) to carry out the necessary expansion, integration, and 
simplification. Here, I present the expansions to order L — 8. 

The choice of expansion parameter affects the compactness 
of the resulting expressions. In the case of Ii, Bessel intro- 



duced a change of variable. 



4e 



(l-e)2' 



Vl + k^- 1 



into Eq. ( l30t to give 

(1 - e)/i(cr) = /" - 2e cos 2cr' + e2 dcr'. 

"'0 



(45) 
(46) 

(47) 



The integrand now exhibits the symmetry that it is invariant 



under the transformation e 



-e and a 



This 



results in a series with half the number of terms (compared 
to a simple expansion in k^). The relation between k and e 
is the same as that between the second eccentricity and third 
flattening of an ellipsoid, e' and n, and frequently formulas for 
ellipsoids are simpler when expressed in terms of n because of 
the symmetry of its definition, Eq. (|2]i. (Bessel undertook the 
task of tabulating the coefficients in the series for Bi{a) for 
some 200 different values of k. This gave him with a strong 
incentive to find a way to halve the amount of work required.) 

The quantity e is 0{f); thus expanding Eq. ( |47l ) to order 
is, asymptotically, equivalent to expanding to order e^. Carry- 
ing out this expansion in e then yields 



A,^il-er\l + 



4*' 



64*' 



256" 



25 ^8 



16384^ 



+•••). 



(48) 
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512" ^ 512" 



" re' + 



16384' 



1280 
2048 



2048" 
4096" 



14336" ^ 

429 8 
r t 



262144' 



(49) 



I use the caesura symbol, ||, to indicate where the series may 
be truncated, at 0{f^), while still giving full accuracy with 
double-precision arithmetic for |/| < 1/150 (this is estab- 
lished in Sect. |9]i. Equation (|49] l is a simple extension of se- 
ries given by iBessell (I1825L §5), except that I have divided 
out the coefficient of the linear term Ai. Bessel's formula- 
tion was used throughout the 19th century an d the serie s given 
here, truncated to order e^, coincides with 'Helmeit' (1880l 
Eq. (5 .5.7)). However, many later works, such as R ainsfordI 
(ll955lEqs. (18)-(19)), use less efficient expansions in k . 
The expansion for I2 proceeds analogously yielding 



^2 = (l-e)(l 
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The expansion of I3 is more difficult because of the pres- 
ence of two parameters / and fc; this presented a problem for 
Bessel — it was a practical impossibility for him to compile 
complete tables of coefficients with two dependencies. Later, 
when the flattening of the earth was known with some preci- 
sion, he might have contemplated compiling tables for a few 
values of /. Instead, Bessel ( 1825, §8) employs a transforma- 
tion to move the dependence on the second parameter into a 
higher order term which he then neglects. The magnitude of 
the neglected term is about 0.000 003" for geodesies stretch- 
ing half way around the WGS84 ellipsoid; this corresponds 
to an error in position of about 0.1 mm. Although this is a 
very small error, there is no need to resort to such trickery 
nowadays because computers can evaluate the coefficients as 
needed. 

Following iHelmertI (1 18801 Eg. (5.8.14)), I expand in n and 
e, both of which are 0{f), to give 
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(53) 



I continue these expansions out to order which is, as noted 
as the end of Sect. ID consistent with expanding the other inte- 
grals to order /^. All the parenthetical terms in Eqs. (|52t and 
(|53] l are functions of n only and so may be evaluated once for 
a given ellipsoid. Note that the coefficient of e' is a terminat- 
ing polynomial of order I in n. This is a curious degeneracy 
of this integral when expressed in terms of n and e. iRainsfordI 
(1 19551 Eqs. (lO)-(l 1)) writes fc^ in Eq. ^ in terms of / and 
cos^ ao and gives a expansion for the integral in powers of /. 
This results in an expansion with more terms. 

The direct geodesic problem requires solving Eq. (l2Ti for a 
in terms of s. (There is a unique solution because ds /da > 0.) 
Equations (|33] | and (|43] |. with j = 1, can be written as 



where 



T = s/{bAl), 



(54) 



(55) 



which shows that finding as a function of s is eq uivalent to 
inverting Eq. (l54l i. This may be accomplished using lLa grange! 
(Em §16) inversion, which gives 



(56) 



where 



(-1)' d'-^B,{a) 



1=1 



l\ 



da'- 



Carrying out these operations with lMaxTm 3 (120091) gives 

oo 

=E^'Si^2ZT, (57) 



where 



C'li 

^12 
Cl3 
Cl4 

C'l, 

^le 
C'lj — 

Cl8 = 



ie_ J.e3 

2^ 32 



16 

29 3 
96 " 



96" ^ 
75 5 
128 



1=1 
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II 516096" 
109167851 
82575360 ' 



2391 6 
2560 " 



1082857 8 
737280 



28223 7 
' 18432" 
733437 
286720 



(58) 



16384" 



iLegendrj ('1806', §13) makes a half-hearted attempt at in- 
verting Eq. (l2Tl i in terms of trigonometric functions of s/h 
(instead of s/{hAi)). Because the period is slightly differ- 
ent from vr, the result is a much more messy expansion than 
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given here. lOrianil disss ') used Lagrange in version to solve 
the distance integral as a series in P. Finally, iHelmerj (Il880i 
Eq. (5.6.8)) carries out the inversion in terms of e (as here) 
including terms to order e^. 

Mo st ot her author s invert Eq. (l2ll iteratively. Both iBessell 
(Il825h and lVincentvl (Il975ah use the scheme given by the first 
line of 

= r-Bi(a«) 

+ A ' 

with ct'"' ~ T. This converges linearly and this is adequate 
to achieve accuracies on the order of 0.1 mm. A superior it- 
erative solution is given by Newton's method which can be 
written as 

Vl + fc^sinVW 

which converges quadratically enabling the solution to full 
machine precision to be found in a few iterations. The sim- 
ple iterative scheme effectively replaces the denominator in 
the fraction above by its mean value Ai. By using the re- 
verted series, I solve for a non-iteratively. This does incur the 
cost of evaluating the coefficients CJ^; however, this cost can 
be amortized if several points along the same geodesies are 
computed. 

I give the expansions for Aj, Cji, and C[i to order 
above. However, these expansions are generated by [Maximal 
(l2009l) with only a few dozen lines of code in about 4 sec- 
onds. The expansions can easily be extended to higher order 
just by changing one parameter in the Maxima code; this has 
been tested by generating the series to order /■^", which takes 
about 13 minutes. Maxima is also used to generate the C-n- 
code for the expansions in GeographicLib. Treating Maxima 
as a preprocessor for the C-n- compiler, the methods presented 
here can be considered to be of arbitrary order In the current 
GeographicLib implementation, the order of the expansion is 
a compile-time constant which can be set to any L < 8. As 
a practical matter, the series truncated at order suffice to 
give close to full accuracy with double-precision arithmetic 
for terrestri al ellipsoids. 

Although lOrianil ( Il806h and lBesseJ d 18251) give expressions 
for general terms in their expansions, most subsequent authors 
are content to work out just a few term s. An exception is 
the work of iLevallois and Dupuvl (Il952h who formulate the 
problem in terms of Wallis integrals where the general term 
in the series is given by a recursion relation, allowin g the se- 
ries to be exten ded to ar bitrary order at run-time (Levalloisl 
[T97O . Chap. 5). iPittmanI (1986) independently derived a sim- 
ilar method. Unfortunately, Pittman uses /3 as the variable of 
integration (instead of a); because the latitude does not change 
monotonically along the geodesic, this choice leads to a loss 
of numerical accuracy and to technical problems in following 
geodesies through vertices (the positions of extrema of the lat- 
itude for the geodesic). 

Because of their w idespread adoption, the expansions given 
by IVuicentvl dl975ah are of special interest. He expresses 



the distance integral in terms of Cn and the longitude inte- 
gral in terms of 1 — ^la. This leads to rather compact series 
when truncated at order However, much of the simplic- 
ity disappears at the next higher order, at which point Vin- 
centy's technique offers no particular advantage. Neverthe- 
less, his procedure does expose the symmetry between k and 
a in Ii{(t) even though his original expansion was in fc^. Be- 
latedly, Vincenty ( 1975a, addendum) discovered the economy 
of Bessel's change of variable, Eq. d45] l. to obtain the same 
expansions for Ai and Cn as given here (truncated at order 
e'^). Incidentally, the key constraint that Vincenty worked un- 
der was that hi s program s should fit onto calculators, such as 
the Wang 720 dVincentvl[T975 b, p. 10), which only had a few 
kilobytes of memory; this precluded the use of higher-order 
expansions and allowed for only a simple (and failure prone) 
iterative solution of the inverse problem. Vincenty cast the 
series in "nested", i.e.. Homer, form, in order to minimize 
program size and register use; I also use the Homer scheme 
for evaluating the series because of its accuracy and speed. 

6. DIRECT PROBLEM 

The direct geodesic problem is to determine <f>2, A12, and 
a2, given 0i, ai and S12; see Fig. [T] The solution starts 
by finding f3i using Eq. dHJ; next solve the spherical triangle 
NEA, to give ao, g\, and uj\, by means of Eqs. (|9l), dT4l l. and 
dTOl ). With ao known, the coefficients Aj, Cji, and C[i may 
be computed from the series in Sect. |5] These polynomials are 
most easily computed using the Horner scheme and the Max- 
ima program accompanying GeographicLib creates the nec- 
essary code. The functions Bj{a) and B[{a), Eq s. d44]i and 
dSTb . may similarly be evaluated for a given a using lClenshawl 
( 119551) summation, wherein the truncated series, 

L 

f{x) = sin^a;, 
1=1 

is computed by determining 

6, = I'' (59) 
I a/ + 2bi+i cos X — &;+2, otherwise, 

and by evaluating the sum as 

f{x) = bi sinx. 

Now si and Ai can be determined using Eqs. d33T l. d35T l. and 
d43b . Compute S2 = si + S12 and find (T2 using Eqs. dSSt and 
d56l ). Solve the spherical triangle NEB to give /32, a2, and uj2 
using Eqs. dTSl ), dT2] i. and dlOl l. Find (/>2 using Eq. ^ again. 
With (72 and W2 given, A2 can be found from Eq. d35l l which 
yields A12 = A2 — Ai. 

Although the reduced length TO12 is not needed to solve the 
direct problem, it is nonetheless a useful quantity to compute. 
It is found using Eqs. dSS, (O, dSll, dllll, (|50|, and dSB- In 
forming /i (cr) — /2 (cr) in Eq. d34l l, I avoid the loss of precision 
in the term proportional to a by writing 

Aia - A2CT = {Ai - l)a - {A2 - l)cr, 
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where Ai — 1, for example, is given, from Eq. (l48T l. by 



^i-l = (l-e)-'(e + 



4" ^64" 



The geodesic scales M12 and Af2i rnay be found similarly, 
starting with Eq. ( l28T l. This completes the solution of the di- 
rect geodesic problem. 

If several points are required along a single geodesic, many 
of the intermediate expressions above may be evaluated just 
once; this includes all the quantities with a subscript "1", and 
the coefficients Aj, Cji, C[i. In this case the determination of 
the points entails just two Clenshaw summations and a little 
spherical trigonometry. If it is only necessary to obtain points 
which are approximately evenly spaced on the geodesic, re- 
place S12 in the specification of the direct problem with the arc 
length on the auxiliary sphere ai2, in which case the conver- 
sion from r to (7 is avoided and only one Clenshaw summation 
is needed. 

For speed and accuracy, I avoid unnecessarily invoking 
trigonometric and inverse trigonometric functions. Thus I 
usually represent an angle 9 by the pair the pair (sin 9, cos 9); 
however, I avoid the loss of accuracy that may ensue when 
computing cos aQ, for example, using — sin^ qq. Instead, 
after finding the sine of ao using Eq. (|9]), I compute its cosine 
with 



cosao = V cos^ a + sin^ a sin^ /?; 
similarly, after finding sin /3 with Eq. (flST l, I use 



i/3 = Vsi 



2 22 
sm ao + cos ao cos a. 



In this way, the angles near the 4 cardinal directions can be 
represented accurately. In order to determine the quadrant of 
angles correctly, I replace Eqs. (fTOt . (fTZt . and (fl4t by 

w = ph(cos(T + i sinofQ sincr), (60) 
a = ph(cos(TCOsao + ^sinao), (61) 
fT = ph(cosacos/3 + isin/3), (62) 

whe re 9 — phfa; + iy) is the phase of a complex num- 
ber dOlveref adl2010[ §1.9(i)), typically given by the library 
function atan2(j/, x). Equations ( |60] l and (|6T| i become inde- 
terminate at the poles, where sin ao — cos ct = 0. However, I 
ensure that uj (and hence A) and a are consistent with their in- 
terpretation for a latitude very close to the pole (i.e., cos is a 
small positive quantity) and that the direction of the geodesic 
in three-dimensional space is correct. In some contexts, the 
solution requires explicit use of the arc length a instead of its 
sine or cosine, for example, in the term Aja in Eq. ( |43] l and 
in determining cr from Eq. (|56l l. 

Geodesies which encircle the earth multiple times can be 
handled by allowing S12 and ai2 to be arbitrarily large. Fur- 
thermore, Eq. ( I6OI 1 allows uj to be followed around the circle in 
synchronism with a; this permits the longitude to be tracked 
continuously along the geodesic so that it increases by +360° 
(resp. —360°) with each circumnavigation of the earth in the 
easterly (resp. westerly) direction. 



The solution for the direct geodesic problem presented here 
is a strai ghtforwa r d exte nsion to higher order of Helmert's 
method ( Helmertl II88OI §5.9), which is largely based on 
iBessell (Il825h . These authors, in common with many more 
recent ones, express the difference of the trigonometric terms 
which arise when the sums, Eq. (l44l i. in Bj{a2) — Bj{ai) are 
expanded as 

sin 2/(72 — sin 2/(71 = 2cos(Z(ct2 + (7i)) sin/(7i2. 

This substitution is needed to prevent errors in the evaluation 
of the terms on the left side of the equation causing large rel- 
ative errors in the difference when using low-precision arith- 
metic. However, the use of double-precision arithmetic ren- 
ders this precaution unnecessary (see Sect.|9); furthermore its 
use interferes with Clenshaw summation and prevents the ef- 
ficient evaluation of many points along a geodesic. 



7. BEHAVIOR NEAR THE ANTIPODAL POINT 

Despite the seeming equivalence of the direct and inverse 
geodesic problems when considered as exercises in ellipsoidal 
trigonometry, the inverse problem is significantly more com- 
plex when transferred to the auxiliary sphere. The included 
angle for the inverse problem on the ellipsoid is A12; however, 
the equivalent angle cji2 on the sphere cannot be immediately 
determined because the relation between A and uj, Eq. ( |23] l, 
depends on th e unknown angle ao- The normal app roach, 
epitomized by iRainsfordI (1 1955b and IVincentvl (Il975ah . is to 
estimate ai and a2, for example, by approximating the el- 
lipsoid by a sphere (i.e., aji2 = A12), obtain a corrected W12 
from Eq. ( |23] | and to iterate until convergence. This procedure 
breaks down if ai and a2 depend very sensitively on 0112, i.e., 
for nearly antipodal points. Before tackling the inverse prob- 
lem, it is therefore useful to examine the behavior of geodesies 
in this case. 

Consider two geodesies starting at A with azimuths ai and 
ai + dai. On a closed surface, they will intersect at some 
distance from A. The first such intersection is the conjugate 
point for the g eodesic and it satisfies mi2 — (for S12 > 0). 
iJacobil d 1 89 l l) showed that the geodesies no longer retain the 
property of being the shortest path beyond the conjugate point 
(iDarbouxl 1 1 894 ?623). For an ellipsoid with small flattening, 
the conjugate point is given by 



<P2 
A2 



Ai 



/vr cos^ 01 cos^ ai + 0(/^), 
■ - /ttcos^i sin^ ai + 0{p). 



The envelope of the geodesies leaving A is given by the lo- 
cus of the conjugate points and, in the case o f an ell i psoid, 
this yields a four-point star called an astroid (iJacobil Il89ll 
Eqs. (16)-(17)), whose equation in cartesian form is, after 
suitable scaling (see below). 



(63) 



whi ch is depicted i n Fig. |4^; see also IJacobil (1189 1'. Fig. 11) 
and lHelmertI d 1 88OI §7.2). The angular extent of the astroid is. 
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30.5 




30.5 



29.5 




179 



179.5 

^2^ 



180 



FIG. 4 Geodesies in antipodal region, (a) The geodesies emanat- 
ing from a point = —30° are shown close in the neighborhood 
of the antipodal point at 02 ~ 30°, A12 = 180°. The azimuths 
Qfi are multiples of 5° between 0° and 180°. The geodesies are 
given in an equidistant cylindrical projection with the scale set for 
02 = 30°. The heavy lines are geodesies which satisfy the shortest 
distance property; the light lines are their continuation. The WGS84 
ellipsoid is used, (b) The geodesic circles on the same scale. The 
heavy (resp. light) lines are for geodesies which have (resp. do not 
have) are property of being shortest paths. The cusps on the circles 
lie on the geodesic envelope in (a). 



to lowest order, 2/7rcos^ 4>i. FigureU^ shows the behavior of 
the geodesic circles in this region. 

Although geodesies are no longer the shortest paths beyond 
the conjugate point (where they intersect a nearby geodesic), 
in general, they loose this property earlier when they first in- 
tersect any geodesic of the same length emanating from the 



same starting point. In the case of the ellipsoid, it is easy to 
establish earlier intersection points. Consider two geodesies 
leaving A with azimuths ai and tt — ai. These intersect at 
|wi2| = TT and 02 = —4>i ^nd, for oblate ellipsoids this in- 
tersection is earlier than the conjugate points. For prolate el- 
lipsoids the corresponding pair of azimuths are ±ai and these 
intersect at |Ai2| — tt, also prior to the conjugate points. Thus, 
an ellipsoidal geodesic is the shortest path if, and only if, 

max(|wi2| , IA12I) < TT. 

The only conjugate points lying on shortest paths are for the 
geodesies with ai — zL^tt for oblate ellipsoids and ai = 
or TT for prolate ellipsoids. Solving the inverse problem with 
end points close to such conjugate pairs presents a challenge 
because tiny changes in end points lead to large changes in the 
geodesic. 

The inverse problem may be solved approximately in the 
case of nearly antipodal points by considering the point B to- 
gether with envelope of geodesies for A (centered at the an- 
tipodes of A); see Fig.|5] The coordinates near the antipodes 
can be rescaled as 



sin(A — Ai — tt) 
AA ' 



y 



sin(/3 + /3i) 



where 



AA = /ttAs cos ,9i , A/3 = cos/3iAA, 



and A3 is evaluated with ao = ^vr — |. In the (x, y) coor- 
dinate system, the conjugate point for the geodesic leaving A 
with ai — zt^TT is at (=Fl, 0) and the scale in the y direction 
compared to a: is 1 + 0( /) . The geodesic through B is tangent 
to the astroid; plane geometry can therefore be used to find its 
direction, using 



y 



cos 9 sin 9 



1 = 0, 



(64) 



where 9 



■^TT ~ OL2~ OL\ 2 

measured anticlockwise from the x axis 



^TT is the angle of the geodesic 
The constant term 
on the left hand side of Eq. (|64] |, 1, reflects the property of 
tangents to the astroid, that the length CD is constant. (The 
astroid is the envelope generated as a ladder slides down a 
wall.) Note that geodesies are directed lines; thus a distinct 
line is given by ^^ — > + tt. The point of tangency of Eq. (|64] i 
with the astroid is 



a:o = — cos 



sm 



which are the parametric equations for the astroid, Eq. 
the line and the astroid are shown in Fig. [5^. The goal now 
is, given x and y (the positio n of B ), to solve Eq. (|64] | for 9. 
I follow the method given bv lVerme ille (2002) for converting 
from geocentric to geodetic coordinates. In Fig.|5ji, COD and 
BED are similar triangles; if the (signed) length BC is k, 
then an equation for k can be found by applying Pythagoras' 
theorem to COD, 

9 2 
(1 + k)2 k2 
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FIG. 5 The inverse geodesic problem for nearly antipodal points, 
(a) The heavy line (labeled 1) shows the shortest geodesic from A 
to B continued until it intersects the antipodal meridian at D. The 
light lines (2-4) show 3 other approximately hemispherical geode- 
sies. The geodesies are all tangent (at their points of conjugacy) to 
the astroid in this figure. The points Q2 and Q3 are the points of 
conjugacy for the geodesies 2 and 3. (b) The solution of the astroid 
equations by similar triangles. 

which can be expanded to give a quartic equation in k, 

+ 2k^ + (1 - - y^)n'^ - 2y'^K -y^ ^0. (65) 

Once K is known, 9 can be determined from the triangle COD 
inFig.lSj), 

e = ph{-x/{l + K)-iy/K). (66) 

The point C in Fig.|5]corresponds to a spherical longitude 
difference of tt. Thus the spherical longitude difference for B, 
UJ12, can be estimated as 

uji2~n + AA; (67) 

1 + K 

compare with lHelmerddlSSOl Eq. (7.3.11)). 

In Appendix |B] I summarize the c losed form solution of 
Eq. (|65ll as given bv lVermeillel(l2002h . I have modified this 



TABLE I The four approximately hemispherical solutions of the 
inverse geodesic problem on the WGS84 ellipsoid for 0i — —30°, 
(j>2 = 29.9°, A12 = 179.8°, ranked by length S12. 



No. 


«i n 


«2(°) 


S12 (m) 


<^i2 n 


mi2 (m) 


1 


161.891 


18.091 


19 989 833 


179.895 


57 277 


2 


30.945 


149.089 


20 010185 


180.116 


24 241 


3 


68.152 


111.990 


20 011887 


180.267 


-22 649 


4 


-81.076 


-99.282 


20 049 364 


180.631 


-68 796 



solution so that it is applicable for all x and y and more stable 
numerically. 

Equation (l65T l has 2 (resp. 4) real roots if B lies outside 
(resp. inside) the astroid. The methods given in Appendix IbI 
(with e set to unity) can be used to determine all these roots, 
K, and Eq. (|66l l then gives the corresponding angles of the 
geodesies at B. All the geodesies obtained in this way are ap- 
proximately hemispherical and that obtained using the largest 
value of K is the shortest path. If B lies on the axes within the 
astroid, then the limiting solutions Eqs. iB6i or ( |B7| l should 
be used to avoid an indeterminate expression. For example, if 
y = 0, substitute the largest k from Eq. (IB7b into Eq. (|66] | to 
give 9 ~ ph(— a; + i-y/max(0, 1 — x^)). 

Figure|5]shows a case where B is within the astroid result- 
ing in 4 hemispherical geodesies which are listed in Table [T] 
ranked by their length. The values given here have been accu- 
rately computed for the case of the WGS84 ellipsoid using the 
method described in Sect. [8] The second and third geodesies 
are eastward (the same sense as the shortest geodesic), while 
the last is westward. As B crosses the boundary of the astroid 
the second and third geodesies approach one another and dis- 
appear (leaving the first and fourth geodesies). Figure |5^ also 
illustrates that the envelope is an evolute of the geodesies; in 
particular. lEisenh^ (1 1 909[ §94) shows that the length of geo- 
desic 2 from A to its conjugate point Q2 exceeds the length of 
geodesic 3 from A to Q3 by the distance along the envelope 
fr om Q3 to Q2; se e also lHelmerll(ll880i §9.2). 

'Schmidt' ('2000') uses Helmert's method for estimating the 
azimuth. Bowring (1996) proposed a solution of the astroid 
problem where he app roximates the 4 arcs of the astroid by 
quarter circles. iRappI (Il993i, Table 1.6, p. 54) gives a simi- 
lar set of hemispherical geodesies to those given in Table [T] 
However, this table contains two misprints: 02 should be 
-40°01'05.759 32" and not -40°00'05.759 32"; for method 
3, ai2 should be 86°20'38.153 06" and not87°20'38.153 06". 

8. INVERSE PROBLEM 

Recall that the inverse geodesic problem is to determine 
S12, «! and oi2 given 0i, 02, and A12. I begin by review- 
ing the solution of the inverse problem assuming that uj\2 is 
given (which is equivalent to seeking the solution of the in- 
verse problem for a sphere). Write the cartesian coordinates 
for the two end points on the auxihary sphere (with unit ra- 
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dius) as A = [cos/3i, 0, sin/3i] and B = [cos/32 cosaji2, 
cos/32sinwi2,sin/32]- A point on the geodesic a small spher- 
ical arc length da from A is at position A + dA where dA — 
[— sin/3i cosai, sinai, cos/3i cos ai] d(T. The azimuth ai 
can be found by demanding that A, B, and dA be coplanar 
or that A X B and A x dA be parallel (where x here denotes 
the vector cross product), and similarly for a2- Likewise, the 
spherical arc length 0-12 is given by ph(A • B + i |A x B|), 
where • denotes the vector dot product. Evaluating these ex- 
pressions gives 

zi = cos l3i sin /32 — sin /3i cos /32 cos 0^12 
+ i cos sina;i2, 
= — sin /3i cos /32 + cos /3i sin /?2 cos W12 
+ i cos f3i sina;i2, 



180 



Z2 



0L\ = phzi, 
a.2 = ph Z2 
CT12 



= Z2 , 

= ph(sin Pi sin /?2 + cos /3i cos P2 cos uji2 + i\zi 



(68) 
(69) 
1). 

(70) 



In order to maintain accuracy when A and B are nearly coin- 
cident or nearly antipodal, I evaluate the real part of zi as 



sin(/32T/3i)± 



sin^ UJ12 sin /3i cos P2 
1 ± cosa;i2 



where the upper and lower signs are for cosa;i2 ^ 0. The 
evaluation of Z2 is handled in the same way. This completes 
the solution of the inverse problem for a sphere. 

In the ellipsoidal case, the inverse problem is just a two- 
dimensional root finding problem. Solve the direct geodesic 
starting at A and adjust ai and S12 (subject to the shortest 
distance constraint), so that the terminal point of the geodesic 
matches B. In order to convert this process into an algorithm, 
a rule needs to be given for adjusting ai and ai2 so that the 
process converges to the true solution. 

The first step in finding such a rule is to convert the two- 
dimensional problem into a one-dimensional root-finding one. 
1 begin by putting the points in a canonical configuration. 



<0, 



< 



< 



< A12 < TT. 



(71) 



This may be accomplished swapping the end points and the 
signs of the coordinates if necessary, and the solution may 
similarly be transformed to apply to the original points. Re- 
ferring to Fig. |3] note that, with these orderings of the coor- 
dinates, all geodesies with ai G [0,7r] intersect latitude (j>2 
with A12 e [0, tt]. Furthermore, the search for solutions can 
be restricted to 02 G [0, ^tt], i.e., when the geodesic ^rif in- 
tersects latitude (j>2- (For 02 = there is a second shortest 
path with a2 G [^tTjTt] if A12 is nearly equal to tt. But this 
geodesic is easily derived from the first.) 

Meridional geodesies are treated as a special case. These 
include the cases A12 = or tt and /?i = —^t^ with any A12. 
This also includes the case where A and B are coincident. In 
these cases, set ai = A12 and 02 = 0. (This value of ai is 
consistent with the prescription for azimuths near a pole given 
at the end of Sect.|5]) 




180 



FIG. 6 The variation of A12 as a function of ai. The latitudes are 
(j>i = -30° and (j>2 = -30°, -25°, -15°, 0°, 15°, 25°, 30°. For 
(fii < and cj>i < 4)2 < — 0i, the curves are strictly increasing, 
while for <^i < and (j>2 — ±<j!>i, the curves are non-decreasing with 
discontinuities in the slopes at qi = 90° (see Fig.|7}. The WGS84 
ellipsoid is used. 



Define now a variant of the direct geodesic problem: given 
4>i, 4)2, subject to Eq. ( fTTT i. and ai, find A12. Proceed as in the 
direct problem up to the solution of the triangle NEA. Find 
(32 from Eq. © and solve the triangle NEB for a2 G [0, ^tt], 
(72, ^2 from Eqs. (|9]l, (fT4l i. and (fTOl l. Finally, determine A12 as 
in the solution to the direct problem. In determining a2 from 
Eq. (|9]l, I use in addition 



+ a/cOs2 tti COS^ Pi + (COS^ ^2 - COS^ ^i) 

cosa2 = — , (72) 

cos P2 

where the parenthetical term under the radical is computed 

by (cos/32 — cos /?! ) (cos /32 + cos/3i) if /3i < — ivr and by 
(sin/3i — sin/32)(sin/3i -I- sin/32) otherwise. It remains to 
determine the value of ai that results in the given value of 

Al2- 

I show the behavior of A12 as a function of ai in Figs. |6]- 
|7] For an oblate ellipsoid and |/32| < — /3i, A12 is a strictly 
increasing function of ai. For P2 — Pi, A12 vanishes for < 
ai < ^tt; for P2 — —Pi, dAi2/dQ!i vanishes for ai = ^'^+- 
Therefore if Pi = P2 = 0, A12 is discontinuous at ai — ^tt 
jumping from to (1 — /)7r. This is the case of equatorial 
end points — if A12 < (1 — /)7r, the geodesic lies along the 
equator, with ai ^ a2 = ^tt. 

Thus with the ordering given by Eq. (iTll . simple root find- 
ing methods, such as binary search or regula falsi, will allow 
ai to be determined for a given A12. Because such meth- 
ods converge slowly, I instead solve for ai using Newton's 
method. 

First, I compute the necessary derivative. Consider a trial 
geodesic with initial azimuth ai. If the azimuth is increased 
to ai + dofi with the length held fixed, then the other end of 
the geodesic moves by mi2 dai in a direction ^tt + q;2. If the 
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FIG. 7 Details of A12 as a function of a\ near the discontinuities 
in the slopes. This shows blow-ups of the two "comer" regions in 
Fig- m (for = —30°). (a) The short line limit, where <f)2 ~ 4>i 
and A12 « 0; here, 02 is in [-30°, -29°] at intervals of 0.1°. At 
4>2 = 4>i = —30° and ai = 90°, the slope changes from zero 
to a finite value, (b) The nearly antipodal limit, where (f>2 ~ —(t)\ 
and A12 « 180°; here, (/>2 is in [29°, 30°] at intervals of 0.1°. At 
4>2 = —4>\ = 30° and ai = 90°, the slope changes from a finite 
value to zero. 



geodesic is extended to intersect the parallel (f)2 once more, the 
point of intersection moves by mi2 dai/ cosa2; see Fig. |8] 
The radius of this parallel is a cos /32, thus the rate of change 
of the longitude difference is 



dAi2 



dai 



^12 



1 



a cos a2 cos P2 



(73) 



where the subscripts on the derivative indicate which quanti- 
ties are held fixed in taking the derivative. The denominator 
can vanish if (32 — \fii \ and a2 — ^tt; in this case, use 



dAi2 



y^^^IL{l^sign{cosa,)), (74) 
sm pi 



for /?2 = ±/?i. ForNewton's method, pick the positive deriva- 
tive, i.e., take 1 =F sign(cosai) = 2, which corresponds to 
ai = 90°± for 02 = ±0i in Fig.|6] 

Newton's method requires a sufficiently accurate starting 
guess for ai to converge. To determine this, I first estimate 
a;i2, the longitude difference on the auxiliary sphere, and then 
find ai using Eq. ( |68] l. To obtain this estimate for a;i2, I dis- 
tinguish three regions; see Fig.|9](the caption gives the precise 
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FIG. 8 Finding dAi2/dQi with 01 and 02 held fixed. 
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FIG. 9 Schematic showing the 5 regions for the inverse problem. 
Region 1 is given by A12 < ^tt and ^2 — Pi < ^tt. Region la is 
given by o"i2 < \/5/ max (0.1, le^l). Re gion 3 is given by cri2 > 
Region 3b is given by j/ > —100(5 and 



7r(l - 3 I/I A3 cos" /3i 
x> -I- 1000^. 



boundaries of the regions). (1) Short lines: If A12 and j32 — Pi 
are reasonably small, then use W12 ~ A12/W1, where wi is 
given by Eq. ( |20l l. (2) Intermediate lines: Assume W12 — A12, 
provided that the resulting 0-12 is sufficiently less than tt. 
(3) Long lines: Analyze the problem using the methods of 
Sect. |7] evaluate a; < and y < 0, and use Eq. Wt\ as an 
estimate of W12. However, if y is very small and —1 < a; < 0, 
then a;i2 is nearly equal to tt and Eq. ( l68T l becomes indeter- 
minate; in this case, estimate ai directly using ai ^ 9 + ^tt, 
with given by Eq. ( |66] l (region 3b in Fig.|9]l. This rule is also 
applied for x slightly less than —1, to ensure that Newton's 
method doesn't get tripped by the discontinuity in the slope of 
Ai2(ai) in Fig.]?]?. 

This provides suitable starting values for ai for use in New- 
ton's method. Carrying out Newton's method can be avoided 
in case 1 above if (T12, Eq. ( iTOb . is sufficiently small (case la in 
Fig- HI), in which case the full solution is given by Eqs. (|68Tl- 
( TTOI ) with S12 = awi<Ji2- This also avoids the problem of 
maintaining accuracy when solving for A12 given 02 ~ 0i 
andai ~ ^tt. Details of the convergence of Newton's method 
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are given in Sect. |9] The boundaries of regions la and 3b in 
Fig. |9] depend on the precision of the floating-point number 
system. This is characterized hy S ~ 1/2^^^ where p is the 
number of bits of precision in the number system and 1 + 6 
is the smallest representable number greater than 1. Typically 
p = 53 and S — 2.2 x 10^^^ for double precision. 

Once a converged value of ai has been found, converged 
values of a2, cr\, and are also known (during the course 
of the final Newton iteration), and Si2 can be found using 
Eq. ( |33] l. The quantities myi, My^, and M21 can also be com- 
puted as in the solution of the direct problem. This completes 
the solution of the inverse problem for an ellipsoid. 

In the following cases, there are multiple solutions to the 
inverse problem and 1 indicate how to find them all given one 
solution. (1) If 01 +01 = (and neither point is at the pole) 
and if ai 7^ a second geodesic is obtained by setting ol\ — 
Q!2 and a2 — oi\. (This occurs when A12 ~ ±180° for oblate 
ellipsoids.) (2) If A12 = ±180° (and neither point is at a pole) 
and if the geodesic in not meridional, a second geodesic is 
obtained by setting a\ — ~a\ and — —ol2- (This occurs 
when the 0i + 02 ~ for prolate ellipsoids.) (3) If A and 
B are at opposite poles, there are infinitely many geodesies 
which can be generated by setting a 1 = ai+7andQ!2 — a^ — 
7 for arbitrary 7. (For spheres, this prescription applies when 
A and B are antipodal.) (4) If S12 = (coincident points), 
there are infinitely many geodesies which can be generated by 
setting ai = ai + 7 and a2 = 0^2 + 7 for arbitrary 7. 

The methods given here can be adapted to return geodesies 
which are not shortest paths provided a suitable starting point 
for Newton's method is given. For example, geodesies 2-\ in 
Table [T] can be found by using the negative square root in the 
equation for cos 02, Eq. (l72T i; and for geodesic 4, solve the 
problem with 27r — A12 = 180.2° as the longitude difference. 
In these cases, the starting points are given by the multiple so- 
lutions of the astroid equation, Eq. (|65] |. Geodesies that wrap 
around the globe multiple times can be handled simil arly. 

Alth ough the published inverse method of Vincentyl 
(Il975a) fails to converge for nearly antipodal points, he did 
give a modificatio n of his method that deals with this case 
(IVincentvl[T975bh . The method requires one minor modifica- 
tion: following his Eq. (10), insert coscr — —\/\ — sin^ g. 
The principal drawback of his method (apart from the lim- 
ited accuracy of his series) is its very slow convergence for 
nearly conjugate points — in some cases, many thousands of it- 
erations are required. In contrast, by using Newton's method, 
the method described here converges in only a few iterations . 
ISodano (1958), starting with the same formulation as given 
here (IBessell 118251: iHelmertL I188O ) derives an approximate 
non-iterative solution for the i nverse problem; however, this 
may fail for antipodal points jRappl 1 19931 §1.3). Sodano's 
justification for his method is illuminating: a non-iterative 
method is better suited to the mechanical and electronic com- 
puters of his day. (See also my comment about Vincenty's 
use of programmable calculators at the end of Sect. |5]) The 
situation now is, of course, completely different: the series of 
Bessel and Helmert are readily implemented on modem com- 
puters and iterative methods are frequently key to efficient and 
accurate computational algorithms. 



TABLE 2 Truncation errors for the main geodesic problems. Ad 
and Ai are approximate upper bounds on the truncation errors for 
the direct and inverse problems. The parameters of the WGS84 
and the SRMmax ellipsoids are used. The SRMmax ellipsoid, a = 
6400 km, / = 1/150, is an ellipsoid with an exaggerated flatten- 
ing introduced by the National Geospatial-lntelligence Agency for 
the purposes of algorithm testing. 
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Ai (m) 
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2.1 X 
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2.1 X 
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-15 


4.1 X 


10- 


-17 


5.2 X 
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-14 
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2.0 X 


10- 
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-25 
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3.0 X 


10- 


-37 


18 


8.5 X 


10- 
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4.2 X 


10~ 
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2.8 X 


10- 
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-48 


1.4 X 


10- 


-53 


1.9 X 


10" 


■41 


2.7 X 


10- 


-47 



9. ERRORS 

Floating-point implementations of the algorithms described 
in Se es. |6] and [8] are included in GeographicLib (Kame^ 
I2OIOI) . These suffer from two sources of error: truncation 
errors because the series in Sect. |5] when truncated at or- 
der L, differ from the exact integrals; and round-off errors 
due to evaluating the series and solving the resulting problem 
in spherical trigonometry using inexact (floating-point) arith- 
metic. In order to assess both types of error, it is useful to be 
able to compute geodesies with arbitrary accuracy. For this 
purpose, I used Maxima's Taylor package to expand the se- 
ries to 30th order and its "bigfloat" package to solve the direct 
problem with 100 decimal digits. The results obtained in this 
way are accurate to at least 50 decimal digits and may be re- 
garded as "exact". 

I first present the truncation errors. I carry out a sequence 
of direct geodesic computations with random 0i, a\, and S12 
(subject to the shortest path constraint) comparing the posi- 
tion of the end point (02 and A12) computed using the series 
truncated to order L with the exact result (i.e., with L — 30), 
in both cases using arithmetic with 100 decimal digits. The 
results are shown in Table |2] which shows the approximate 
maximum truncation error as a function of L < 20. The quan- 
tity Ad gives the truncation error for the method as given in 
Sect.|6l this scales as On the other hand, Ai is the trunca- 
tion error where, instead of solving a in terms of s using the 
truncated reverted series, Eq. (ISFt , I invert the truncated se- 
ries for s, Eqs. ( l33T l and ( |43] |. to give a in terms of s. (This is 
done "exactly", i.e., using Newton's method and demanding 
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convergence to 100 decimal places.) This is representative of 
the truncation error in the solution of inverse geodesic prob- 
lem, because the inverse problem does not involve determin- 
ing CT in terms of s. Aj scales approximately as {^f)^- For 
comparison, the truncation errors for Vincenty's algorithm are 
9.1 X IQ-s ni and 1.5 x lO^^ m for the WGS84 and SRMmax 
ellipsoids; these errors are about 2.5 times larger than Ad for 
L = 3 (the order of Vincenty's series). 

1 turn now to the measurement of the round-off errors. 
The limiting accuracy, assuming that the fraction of the 
floating-point representation contains p = 53 bits, is about 
20 000 km/2^^ « 2 nm (where 20 000 km is approximately 
half the circumference of the earth). From Table |2] the choice 
L = 6 ensures that the truncation error is negligible compared 
to the round-off error even for / = 1/150. 1 assembled a large 
set of exact geodesies for the WGS84 ellipsoid to serve as test 
data. These were obtained by solving the direct problem using 
Maxima using the protocol described at the beginning of this 
section. All the test data satisfies 0-12 < 180°, so that they 
are all shortest paths. Each test geodesic gives accurate values 
for (pi, ai, 4>2, A12, S12, (J12, and mi2. In this list <j)i, ai, 
and S12 are "input" values for the direct problem. The other 
values are computed and then rounded to the nearest 0.1 pm 
in the case of mi2 and (10^^^)° in the case of the angles. The 
test data includes randomly chosen geodesies together with 
a large number of geodesies chosen to uncover potential nu- 
merical problems. These include nearly meridional geodesies, 
nearly equatorial ones, geodesies with one or both end points 
close to a pole, and nearly antipodal geodesies. 

For each test geodesic, 1 use the floating-point implementa- 
tions of the algorithms to solve the direct problem from each 
end point and to solve the inverse problem. Denoting the re- 
sults of the computations with an asterisk, 1 define the error in 
a computed quantity xhy 5x = x* —x. For the direct problem 
computed starting at the first end point, 1 compute the error in 
the computed position of the second end point as 

|p2<502 + i cos 4>2l'2SX2 I ■ 

I convert the error in the azimuth into a distance via 

a \Sa2 ~ (5Ai2 sin 02 1 , 

which is proportional to the error in the direction of the geo- 
desic at B in three dimensions and accounts for the coupling 
of a2 and A12 near the poles. 1 compute the corresponding 
errors when solving the direct problem starting at the second 
end point. 

For the inverse problem, 1 record the error in the length 

|<^S12|. 

I convert the errors in the azimuths into a length using 

max(|(5Q!i| , |fo2|) |mi2| ; 

the multiplication by the reduced length accounts for the sen- 
sitivity of the azimuths to the positions of the end points. An 
obvious example of such sensitivity is when the two points 
are close to opposite poles or when they are very close to each 



TABLE 3 Two close geodesies. The parameters of the WGS84 el- 
lipsoid are used. 



Case 1 Case 2 

01 -30° 

02 30° (30 - 

4 X 10-^^)° 
A12 179.477 019 999 756 66° 

ai 90.000 008° 90.001 489° 

02 89.999 992° 89.998 511° 
si2 19 978 693.309 037 086 m 
cri2 180° 180.000 000° 

mi2 1.1 nm 51 /im 



other A less obvious case is illustrated in Table |3] The second 
end points in cases 1 and 2 are only 0.4 nm apart; and yet the 
azimuths in the two cases differ by —5.3" and the two geo- 
desies are separated by about 160 m at their midpoints. (This 
"unstable" case finding the conjugate point for a geodesic by 
solving the direct problem with a\ — 90° and spherical arc 
length of (T12 = 180°.) 

This provides a suitable measure of the accuracy of the 
computed azimuths for the inverse problem. 1 also check their 
consistency. If S12 > a, I determine the midpoint of the com- 
puted geodesic by separate direct geodesic calculations start- 
ing at either end point with the respective computed azimuths 
and geodesic lengths ±^Si2 and I measure the distance be- 
tween the computed midpoints. Similarly for shorter geode- 
sies, s\2 < a, 1 compare the computed positions of a point 
on the geodesic a distance a beyond the second point with 
separate direct calculations from the two endpoints and repeat 
such a comparison for a point a distance a before the first end 
point. 

In this way, all the various measures of the accuracy of 
the direct and inverse geodesic are converted into comparable 
ground distances, and the maximum of these measures over a 
large number of test geodesies is a good estimate of the com- 
bined truncation and round-off errors for the geodesic calcu- 
lations. The maximum error using double precision (p — 53) 
is about 15 nm. With extended precision (p — 64), the error 
is about 7 pni, consistent with 1 1 bits of additional precision. 
The test data consists of geodesies which are shortest paths, 
the longest of which is about 20 000 km. If the direct problem 
is solved for longer geodesies (which are not therefore short- 
est paths), the error grows linearly with length. For example, 
the error in a geodesic of length 200 000 km that completely 
encircles the earth 5 times is about 150 nm (for double preci- 
sion). 

Another important goal for the test set was to check the con- 
vergence of the inverse solution. Usually a practical conver- 
gence criterion for Newton's method is that the relative change 
in the solution is less than 0(\/J); because of the quadratic 
convergence of the method, this ensures that the error in the 
solution is less that 0{8). However, this reasoning breaks 
down for the inverse geodesic problem because the deriva- 
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TABLE 4 Ellipsoidal trigonometry problems. Here, ^, (, and 9 are 
the three given quantities. The "notes" column gives the number as- 
signed bV|Oriani ( 1810 , p. 48) in his "index of spheroidal problems" 
and by iPuissan 521) in his enumeration of solutions. See 

the text for an explanation of the other columns. 



No. 






Ref. 




5,C 


Notes 


1 


CD^ (Ti (J)o 










1, 1 


2 


01, O?!, 02 










3, 12 


3 


01, Ql, Si2 


0"12 




Eq. ^ 


5,2 


4 


01, ai, Ai2 


0"12 




Eq. (|76l( 


11,6 


5 


01, 02, Sl2 


Ql 


1 


Eq. ((TTj 


7,3 


6 


01, 02, Ai2 


ai 


1 


Eq. ^ 


13,7 


7 


01, Q2, Si2 


Ql 


2 


Eq. (I?! 


8,4 


8 


01, 02, Ai2 


Ql 


2 


Eq.d 


loll 


14, 8 


9 


01, Sl2, Al2 


Ql 


3 


Eq. (IB 


19, 11 


10 


ai,si2, Ai2 


01 


3 


Eq.d 


n 


17, 10 


11 


ai, Q!2, Sl2 


01 


2 


Eq. ([83 


10,5 


12 


ai, a2, Ai2 


01 


2 


Eq. ^ 


16,9 



live of Ai2 with respect to ai can become arbitrarily small; 
therefore a more conservative convergence criterion is used. 
Typically 2-4 iterations of Newton's method are required. A 
small fraction of geodesies, those with nearly conjugate end 
points, require up to 16 iterations. No convergence failures 
are observed. 



10. ELLIPSOIDAL TRIGONOMETRY 

The direct and inverse geodesic problems are two examples 
of solving the ellipsoidal triangle NAB in Fig. [T] given two 
sides and the included angle. The sides of this triangle are 
given by NA = aE{^Tr - ^i,e), NB = aE{^n - /?2,e), 
and AB = si2 and its angles are NAB = a\, NBA = tt — 
(32, and ANB = Ai2. The triangle is fully solved if 0i, ol\, 
02, «2, Si2, and Ai2 are all known. The typical problem in 
ellipsoidal trigonometry is to solve the triangle if just three of 
these quantities are specified. Considering that A and B are 
interchangeable, there are 12 distinct such problems which are 
laid out in Table H) (In plane geometry, there are four distinct 
triangle problems. On a sphere, the constraint on the sum 
of the angles of a triangle is relaxed, leading to six tiiangle 
pr oblems .) 

lOrianil (llSlOi p. 48) and lPuissanj (Il83ll p. 521) both gave 
similar catalogs of ellipsoidal problems as Table |4] Here (and 
in Sect. [TTT i. I do not give the full solution of the ellipsoidal 
problems nor do I consider how to distinguish the cases where 
there may be 0, 1, or 2 solutions. Instead, I indicate how, in 
each case, an accurate solution may be obtained using New- 
ton's method assuming that a sufficiently accurate starting 
guess has been found. This might be obtained by approximat- 
ing the elhpsoid by a sphere and using spherical trigonome- 
try dTodhunteil 1 1 87 1[ Chap. 6) or by using approximate ellip- 



soidal methods ('RappVl99l', §6). In Sect. [131 1 also show how 
the gnomonic projection may be used to solve several ellip- 
soidal problems using plane geometry. 

In treating these problems, recall that the relation between 
(\) and /? is given by Eq. (O and depends only on the eccentric- 
ity of the ellipsoid. On the other hand, the relations between 
si2 and Ai2 and the corresponding variables on the auxiliary 
sphere, (j\2 and wi2, depend on the geodesic (specifically on 
ao)- 

In Table m ^, C, ^ are the given quantities. In problems 1 
and 2, the given quantities are all directly related to corre- 
sponding quantities for the tiiangle on the auxiliary sphere. 
This allows the auxiliary triangle to be solved and the ellip- 
soidal quantities can then be obtained. (Problem I was used 
in solving the inverse problem in Sect. |8]) 

Problem 3 is the direct problem whose solution is given 
in Sect. |6] Here, however, I give the solution by Newton's 
method to put this problem on the same footing as the other 
problems. The solution consists of treating u\2 (the col- 
umn labeled t/j) as a "control variable". Assume a value for 
this quantity, solve the problem with given 0i,q:i,<ti2 (i.e., 
^, C, -0) for si2 (i.e., ff) using Eq. ( [33] l. The value thus found 
for si2 will, of course, differ from the given value and a better 
approximation for a\2 is found using Newton's method with 



dsi2 



dCTl2 



(75) 



The equation for the derivative needed for Newton's method is 
given in the column labeled d6'/d?/'|^ ^ in the table. Problem 



4 is handled similarly except that Eq. ( f35] l is used to give A12 
and the necessary derivative for Newton's method is 



dAi 



dCTl2 



U)2 Sllia2 
COS {^2 



(76) 



Problems \-4 are the simplest ellipsoidal trigonometry 
problems with and a specified at the same point, so that it is 
possible to determine olq which fixes the relation between S12 
and (T12 and between A12 and bj\2. In the remaining problems 
it is necessary to assume a value for 0i or oi\ thereby reducing 
the problem to one of the reference problems 1-3 (the column 
labeled "Ref."). Thus, given C, assume a value -0^°-' for -0; 
solve the reference problem ^, V''^*-' to determine 6''^*); find a 
more accurate approximation to ?/; using 



d6'(^) 



and iterate until convergence. The remaining derivatives are 

dsi2 



dai 
dAi2 



dai 
dsi2 



01,02 



dai 



mi2 tan 02, 

■mi_2 1 
a cos ai cos ^2 ' 

mi2 tan a2 



aw2 



tan a.\ tan ^2 cos a2 



(77) 
(78) 
(79) 



16 



dA 



12 



dAi2 



dai 
dAi2 



f>l ,C(2 



01,S12 



il,Sl2 



myija W2tana2 
cos Q!2 cos /32 tan ai sin fii 

mi2 cos a2 
a cos /32 ' 

3 

(1 - /) sinai 

/ cosai cosq;2 



(80) 
(81) 



dsi2 



^ cos P- 



= a 



ai ,02 



1-/' 



dA 



42 



ai ,a2 



COS /?2 

A'^12 tan a2 
sin ai 

tan Pi 

H ^-^^ 

cos a2 tan /32 wi 

cos ai 



12 



W2 



(82) 



(83) 



Ni2 sec a2 
^sinaicos/?! sinaiCOs/32 
^ tan /3i tan a2 W2 \ 



sin/32 



(84) 



where 



7V- 



12 



Mi2 



(mi2 /a) cos ai tan Pi 

Wl 



The choice of ip in these solutions is somewhat arbitrary; 
other choices may be preferable in some cases. These formu- 
las for the derivatives are obtained with constructions similar 
to Fig.|8] Equations (I82]i-(l84li involve partial derivatives taken 
with ai held constant; the role of N12 in these equation can 
be contrasted with that of M12 as follows. Consider a geo- 
desic from A to B with length S12 and initial azimuth ai. 
Construct a second geodesic of the same length from A' to 
B' where A' is given by moving a small distance dt from A 
in a direction ai + ^tt. If the initial direction of the second 
geodesic is a'l — ai (resp. parallel to the first geodesic), then 
the distance from B to B' is N12 dt (resp. M12 dt). (Because 
meridians converge, two neighboring geodesies with the same 
azimuth are not, in general, parallel.) 

Problem 6 is the geodesic inverse problem solved in Sect. |8] 
and 1 have repeated Eq. ( l73T l as Eq. (iTSTl. Problem 7 is the 
"retro-azimuthal" problem for which iHinksl (1 19291 ) gives an 
interesting application. For many years a radio at Rugby trans- 
mitted a long wavelength time signal. Hinks' retro-azimuthal 
problem is to determine the position of an unknown point with 
knowledge of the distance and bearing to Rugby. 



11. TRIANGULATION FROM A BASELINE 

The ellipsoidal triangle considered in Sect. [TO]is special in 
that one of its vertices is a pole so that two of its sides are 
meridians. The next class of ellipsoidal problems treated is 
solving a triangle ABC with a known baseline, see Fig.fTOh. 
Here, the goal is to determine the position of C if the positions 
if A and B are known, i.e., if 0i, 02, and A12 are given (and 
hence, from the solution of the inverse problem for AB, the 
quantities ai, a2, and S12 are known). This corresponds to a 



(a) C 




(b) C 



FIG. 10 Two ellipsoidal triangle problems: (a) triangulating from a 
baseline and (b) rectangular geodesic (or oblique Cassini-Soldner) 
coordinates. 



TABLE 5 Ellipsoidal triangulation problems. In these problems, 
the positions of A and B (i.e., 0i, <j>2, and A12) are given and the 
position of C is sought. The quantities ^ and are the additional 
given quantities. 



No. 






dC 

dtp 


e 





Q!l(3), S13 








I 


S13, S23 


ai(3) 


Eq. O 




2 


«1(3), «2(3) 


Sl3 


Eq. O 


Hi 


3 


Ql(3)i S23 


Sl3 


Eq. lETll 


4 


Sl3,73 


"1(3) 


Eq. O 


3 


5 


ai(3),73 


Sl3 


Eq. 


9} 



class of triangulation problems encountered in field surveying: 
a line AB is measured and is used as the base of a triangula- 
tion network. However, in the present context, A need not be 
visible from B and 1 consider both the problems of triangula- 
tion and trilateration. (In addition, remember that the angles 
measured by a theodolite are not the angles between geodesies 
but between normal sections.) 

Because the problems entail consideration of more than a 
single general geodesic, it is necessary to generalize the nota- 
tion for azimuthal angles to make clear which geodesic line is 
being measured. I define cti(j) as the azimuth of the geodesic 
line passing through point i where j is some other point on 
the same line. Each geodesic line is assigned a unique direc- 
tion, indicated by arrows in the figures, and all azimuths are 
forward azimuths (as before). 

With one side specified and two additional quantities 
needed to solve the triangle, there are 10 possible problems to 
solve. Of these, six are distinct (considering the interchange- 
ability of A and B) and are listed in Table |5] For the angle 
at C, I assume that 73 — '^3(1) — ^3(2)? the difference in the 
bearings of A and B, is given; this is included angle at C in 
Fig. [TOb . In other words, 1 do not presume that the direction 
of due north is known, a priori, at C. This is the common 
situation with theodolite readings and it also includes the im- 
portant case where 7 is required to be i^Tr which allows the 



17 



shortest distance from a point to a geodesic to be determined. 

Problem is just the direct geodesic problem (problem 3 of 
Sect. [Toll. The remaining 5 problems may be solved by New- 
ton's method, in a similar fashion as in Sect. (TO] as follows: 
replace the second given quantity ( by one of the unknowns ip; 
estimate a value of ip; solve the problem with ^ and tp (which, 
in each case, is a direct geodesic problem from A) to deter- 
mine a trial position for C; solve the inverse geodesic problem 
between B and the trial position for C to obtain a trial value 
for (; update the value of ip using Newton's method so that the 
resulting value of ( matches the given value. The derivatives 
necessary for Newton's method are given in Eqs. 



ds 



23 



dai 



(3) 



!>l,02,Al2,Sl3 



2(3) 



dsi3 

ds23 



-TOi3Sm73, 

1 . 



dsi3 
d73 



^1,02,^12,0:1(3) 



01,02,Ai2,Q;i(3) 



m23 

COS 73, 



■sm73, 



(85) 
(86) 
(87) 



dai(3) 
d73 



ds 
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!>l,02,Al2,Sl3 



i'l,02,Ai2,ai(3) 



M3i-M32^COS73, (88) 
"^23 



M32 

sin 73. 

'7^23 



(89) 



The triangulation problems 1, 2, and 3, entail specification 
of either the distance or the bearing to C from each of A and 
B. These can also be solved by generalizing the method given 
bv lSiober3(l2002l §5) for the solution of problem 2. The tech- 
nique is to treat the latitude of C, (/13, as the control variable. 
Thus, start with an estimate for (j)^; if the bearing of (resp. dis- 
tance to) C from a base point is given, then solve the inter- 
mediate problem 4^, 03, ans) (resp. 5,^3) where i = 1 or 2 
for base points A or B, i.e., problem 1 (resp. 5) in Sect. [TOlto 
give Ai3; and evaluate A12 = A13 ~ A23. Now adjust ^3 us- 
ing Newton's method so that the A12 matches the known value 
using 



dA, 



i3 



tan 



"3(1 



1(P3 

dAi3 



1-/ cos/ 



Wo 



cot a3(j 



1 - / cos /33 



(90) 
(91) 



This method of solution essentially factors the problem into 
two simpler problems of the type investigated in Sect. [10] 

A similar approach can also be applied to problem 5. Guess 
a value of S13, solve problems 3 and 1 of Sect. [TOjto deter- 
mine successively the positions of C and B. Adjust S13 using 
Newton's method so that the correct value of A12 is obtained, 
which requires the use of the derivative 



dAi 



dsi3 



P2,ai(3),73 



M32 sin 73 sec a2 (3) 
a cos /32 



(92) 



Knowledge of reduced length and the geodesic scale allow 
errors to be propagated through a calculation. For example. 



if the measurements of q;i(3) and 02(3) are subject to an in- 
strumental error Sa, then the error ellipse in the position of 
C when solving problem 2 will have a covariance which de- 
pends on 77113 ^oi^ "7,23 <5a, and 73. The effects of errors in the 
positions of A and B and in the measurements of S13 and S23 
c an be similarly estimated. 

iRNAVl (I2OO7I §§ A2.4.1-3) also presents solution for prob- 
lems 1-3. However, these use the secant method and so con- 
verge more slowly that the methods given here. 



12. GEODESIC PROJECTIONS 

Several map projections are defined in terms of geod esies. 
In the azimuthal equidistant projection dSnvdeil Il987i §25) 
the distance and bearin g from a centr al point A to an arbitrary 
point B is preserved. GaussI (Il902[ §19) lays out the prob- 
lem for a general surface: the point B is projected to plane 
cartesian coordinates. 



Si2smai, 7/ = S12 cos Qfi. 



iBagratunil ( 1967[ §1 6) calls these "geodetic polar coordi- 



nates". Gauss IT9O2I §15) proves that the geodesies (lines of 
constant ai) and the geodesic circles (lines of constant S12), 
which, by construction, intersect at right angles in the projec- 
tion, also intersect at right angles on the ellipsoid (see Sect.[3]|. 
The scale in the radial direction is unity, while the scale in 
the azimuthal direction is 512/7)712; the projection is confor- 
mal only at the origin. The forward and reverse projections 
are given by solving the inverse and direct geodesic problems. 
The entire ellipsoid maps to an approximately elliptical area, 
with the azimuthal scale becoming infinite at the two bound- 
ary points on the x axis. The projection can be continued 
beyond the boundary giving geodesies which are no longer 
shortest lines and negative azimuthal scales. ISnvded (1 19871 
p. 197) gives the formulas for this projection for the ellipsoid 
only for the case where the center point is a pole. For example, 
if A is at the north pole then the projection becomes 

S12 ^ aE{^7r ~ I32,e), 77112 = a cos /32. 

However, the method given here is applicable for any center 
point. The projection is useful for showing distances and di- 
re ctions f r om a c entral transportation hub. 

iGaussI (I19O2I §23) also describes another basic geode- 
sic pro jection, called "right-angle spheroidal coordinates" by 
[Bagrat uni ( 1967, § 17). Consider a reference geodesic passing 
through the point A at azimuth ai(3). The reverse projec- 
tion, B, of the point x, y is given by the following operations 
which are illustrated in Fig.[Tob: resolve the coordinates into 
the directions normal and parallel to the initial heading of the 
reference geodesic, 

x' — cosai(3)X — sinai(3-)?/, 
y' = sin ai(^3)X + cos ai(3)?/; 

starting at A proceed along the reference geodesic a dis- 
tance y' to C; then proceed along the geodesic with azimuth 
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Q^3(2) = 0^3(1) + ^TT a distance x' to point B. (Here, the "for- 
ward" direction on the geodesic CB is to the right of the ref- 
erence geodesi c AC which is the opposite of the convention 
in Sect. [lU) Gauss ('1902', §16) proves that geodesies (lines 
of constant y') and the geodesic "parallels" (lines of constant 
x') intersect at right angles on the ellipsoid. At B, the scale 
in the x' direction is unity while the scale in the y' direction 
is l/A/32 which is unity on x' = 0; thus the projection is con- 
formal on x' — 0. The definition of the mapping given here 
provides the prescription for carrying out the reverse proj ec- 
tions. The forward projection is solved as follows: determine 
the point C on the reference geodesic which is closest to B 
(problem 5 of Sect. II lb: set x' to the distance CB signed pos- 
itive or negative according to whether B is to the right or left 
of the reference geodesic; set y' to the distance AC signed 
positive or negative according to whether C is ahead or behind 
A on the reference geodesic; finally transform the coordinate 
frame 



x = cosair.^ia;' + sinaijsjy' 



V 



sinQ!i(3)x' + cosai(3)?/' 



In the case where the reference geodesic is the equator, 
the projection is the ellipsoidal generalizat ion of t he so- 
called "equidistant cylindrical" projection (ISnvderl Il987[ 
§12). Solving for the point C is trivial; the coordinates are 
given by 

x = aAi2, y a(£;(e)-£;(i7r-/?2,e)), Af32 = cos/32, 

where -E'(fc) i s the corn p lete e lhptic integral of the 
second kind jOlver ef a/.l I2OIOI §19.2(ii)); see also 
iBugavevskiv and Snvdei ( 1995 , §2.1.4). The whole ellipsoid 
is mapped to a rectangular region with the poles mapped to 
lines (where the scale in the x direction is infinite). If the refer- 
ence geodesic is a "central" merid ian, the projection is called 
"Cassini-Soldner' ' (ISnvdeil[T987l §13) and C may most sim- 
ply be found by finding the midpoint of the geodesic BD 
where D is the reflection of B in the plane of the central 
meridian. This allows the Cassini-Soldner mapping to be be 
solved accurately for the whole eUipsoid , in contrast to the se- 
ries method presented in ISnvderl (Il987t p. 95) which is only 
valid near the central meridian. The ellipsoid maps to an ap- 
proximately rectangular region with the scale in the y direc- 
tion divergent where the equator intersects the boundary. The 
Cassini-Soldner was widely used for large-scale maps until 
the middle of the 20th century when it was alm ost entirely re- 
placed by the transverse Mercator projection (I Snvderl Il987l 
§8). Tasks such as navigation and artillery aiming were much 
more easily accomplished with a conformal projection, such 
as transverse Mercator, compared to Cassini-Soldner with its 
unequal scales. The general case of this mapping may be 
termed the "oblique Cassini-Soldner" projection. Because the 
reference geodesic is not closed in this case, it is not conve- 
nient to use this projection for mapping the entke ellipsoid 
because there may be multiple candidates for C, the position 
on the reference geodesic closest to B. 

The doubly equidistant projection has been used in small- 
scale maps to minimize the distortions of large land masses 




FIG. 1 1 The doubly equidistant projection. Also shown are portions 
of the geodesic ellipse and hyperbola through C. 



dBugavevs kiv and Snvdeii Il995l §7.9). In this projection, the 
distances to an arbitrary point C from two judiciously cho- 
sen reference points A and B are preserved; see Fig. [TT] The 
formulas are usually given for a sphere; however, the general- 
ization to an ellipsoid is straightforward. For the forward pro- 
jection, fix the baseline AB with A and B separated by syi, 
solve the inverse geodesic problems for AC and BC and use 
the distances §13 and §23 together with elementary trigonome- 
try to determine the position of C in the projected space; there 
are two solutions for the position of C either side of the base- 
line; the desired solution is the one that lies on the same side 
of the baseline as C on the ellipsoid. The reverse projection is 
similar, except that the position of C on the ellipsoid is deter- 
mined by solving problem 1 in Sect.[TT] The projection is only 
well defined if 71 and 72 are the same sign (consistent with a 
planar triangle). This is always the case for a sphere; the en- 
tire sphere projects onto an ellipse. However, on an ellipsoid, 
the geodesic connecting A and B is not closed in general and 
thus does not divide the ellipsoid into two halves. As a conse- 
quence, there may be a portion of the ellipsoid which is on one 
side of the baseline geodesic as seen from A but on the other 
side of it as seen from B\ such points cannot be projected. For 
example if A = (35°N, 40°E) and B = (35°N, 140°E), then 
the point (43.5°S, 60.5°W) cannot be projected because it is 
north of the baseline as seen by A but south of the baseline 
relative to B 

The scales of the doubly equidistant projection can be de- 
termined as follows. By construction, geodesic ellipses and 
hyperbolae, defined by u± — 5(313 ± S23) — const., map 
to ellipses a r id hyp erbolae under this projection; see Fig. [TTI 
IWeingartenI (1186 3') establi shes these re sults for geodesic el- 
lipses and hyperbolae (Eis enhartlll909[ §90): they are orthog- 
onal; the geodesic hyperbola through C bisects the angle 73 
made by the two geodesies from A and B\ and the scales in 
the u± directions are cos ^73 and sin ^73, respectively. The 
same relations hold, of course, for the projected ellipses and 
hyperbolae, except that the angle ACB takes on a different 
(smaller) value 73. Thus, the elliptic and hyperbolic scale fac- 



19 



c 




c 



A B 



FIG. 12 The construction of the spheroidal gnomonic projection as 
the limit of a doubly azimuthal projection. 

tors for the double equidistant projection may be written as 

cos sin 
cos i7^ ' sin ^7^ ' 

respectively. Evaluating 73 using the cosine rule for the plane 
triangle ABC gives 

2^/513523 cos ^73 2^513523 sin ^73 

V(S13 + S23)^ - sj^ \Js\2 - (Sl3 - S23)^ 

The results generalize those of lCoxl (1 19461 Il95lh for the pro- 
jection of a sphere. In the limit S\2 — > 0, this projection re- 
duces to the azimuthal equidistant projection and these scales 
reduce to the radial scale, 1, and the azimuthal scale, S13 /mi3. 

13. SPHEROIDAL GNOMONIC PROJECTION 

The gnomonic projection of the sphere, which is obtained 
by a central projection of the surface of the sphere onto a tan- 
gent plane, has the pr operty that a ll geodesies on the sphere 
map to straight lines (ISnvderi Il987i, §22). Such a projection 
is impossible f or an ellipsoid b ecause it does not have con- 
stant curvature (lBeltramilll865l) . However, a spheroidal gen- 
eralization of the gnomonic projections can be constructed for 
which geodesies are very ne arly straight. First recall that the 
doubly azimuthal projection ( Bugavev skiv and Snvdeil 1 19951 
§7.8) of the sphere, where the bearings from two points A and 
B to C are preserved, gives the gnomonic projection which is 
compressed in the direction parallel to AB. In the limit as B 
approaches A, the pure gnomonic projection is recovered. 

The construction of the spheroidal gnomonic projection 
proceeds in the same way; see Fig. [12] Draw a geodesic 
BC such that it is parallel to the geodesic AC at B. Its ini- 
tial separation from AC is Si2sin7i; at C", the point clos- 
est to C, the separation becomes M13S12 sin 71 (in the limit 



S12 0). Thus the difference in the azimuths of the geode- 
sies BC and BC at B is (Afi3/mi3)si2 sin 71, which gives 
71 + 72 = — (Mi3/rni3)si2 sin 71. Now, solving the planar 
triangle problem with 71 and 72 as the two base angles gives 
the distance AC in the projected space as 77113/1/13. 

Thus leads to the following specification for the spheroidal 
gnomonic projection. Let the center point be A; for an arbi- 
trary point B, solve the inverse geodesic problem between A 
and B\ then point B projects to the point 

a; = psinai, y = /9C0sai, p — myilM^vi^ (93) 

the projection is undefined if A/12 < 0. In the spherical limit, 
this becomes th e standard gnomonic projection, p ^ a tan (T12 
(ISnvdeii[T98 7', p. 165). The azimuthal scale is l/A/12 and the 
radial scale, found by computing <:\p/dsi2 and using Eq. ( |29] |, 
is l/A/12; the projection is therefore conformal at the origin. 
The reverse projection is found by ai — ph(?/ + ix) and 
by solving for S12 using Newton's method with d/9/dsi2 — 
l/A/^2 the radial scale). Clearly the projection preserves 
the bearings from the center point and all lines through the 
center point are geodesies. Consider now a straight line BC 
in the projection and project this line on the spheroid. The 
distance that this deviates from a geodesic is, to lowest order, 

h=^^{VK-t)t, (94) 

where I is the length of the geodesic, K is the Gaussian curva- 
ture, and t is the perpendicular vector from the center of pro- 
jection to the geodesic. 1 obtained this result semi-empirically: 
numerically, I determined that the maximum deviation was 
for east-west geodesies; 1 then found, by Taylor expansion, 
the deviation for the simple case in which the end points are 
equally distant from the center point at bearings zta; finally, 
I generalized the resulting expression and confirmed this nu- 
merically. The deviation in the azimuths at the end points is 
about 4/i/^ and the length is greater than the geodesic dis- 
tance by about ^h? /I. For an ellipsoid, the curvature is given 
by Eq. ( |25] ), which gives 

VK ^ -^6^(1 - e^sm^(f>f/^cos(l)sm(f>; (95) 

the direction of V/C is along the meridian towards the equator. 
Bounding h over all the geodesies whose end-points lie within 
a distance r of the center of projection, gives (in the Umit that 
e and r are small) 

\h\<l^^r. (96) 

The limiting value is attained when the center of projection 
is at (p — ±45° and the geodesic is running in an east-west 
direction with the end points at bearings ±45° or ±135° from 
th e center. 

iBowrind (1 19971) and lWilliarnsI (Il997h have proposed an al- 
ternate ellipsoidal generaUzation of the gnomonic projection 
as a central projection of the ellipsoid onto a tangent plane. In 
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FIG. 13 The coast line of Europe and North Africa in the ellipsoidal 
gnomonic projection with center at (45°N, 12°E) near Venice. The 
graticule lines are shown at multiples of 10°. The two circles are 
centered on the projection center with (geodesic) radii of 1000 km 
and 2000 km. The data for the coast lines is taken from GMT 
fWessel and Smith.,201Q) at "low" resolution. 



This can be solved using the ellipsoidal gnomonic projection 
as follows. Guess an intersection point O^^^^ and use this as 
the center of the gnomonic projection; define a, b, c, d as the 
positions of A, B, C, D in the gnomonic projection; find the 
intersection of of AB and CD in this projection, i.e., 

(z • c X d)(b — a) — (z • a X b)(d — c) 
z • (b — a) X (d — c) 

where ' indicates a unit vector (a = a/a) and z = x x y is 
in the direction perpendicular to the projection plane. Project 
o back to geographic coordinates O'^^ and use this as a new 
center of projection; iterate this process until O^'^ — O*^'^^^ 
which is then the desired intersection point. This algorithm 
converges to the exact intersection point because the mapping 
projects all geodesies through the center point into straight 
lines. The convergence is rapid because projected geodesies 
which pass near the center point are very nearly straight. Prob- 
lem 5 of Sect. [TT] can be solving using the gnomonic projec- 
tion in a similar manner If the point O on AB which closest 
to C is to be found, the problem in the gnomonic space be- 
comes 

c • (b — a)(b — a) — (z • a X b)z x (b — a) 



such a mapping, great ellipses project to straight lines. Em- 
pirically, I find that the deviation between straight lines in this 
mapping and geodesies is 

\h\ < |-r. 
2 a 



iLetovaF tsevi d 1 9631) suggested another gnomonic projection in 
which normal sections through the center point map to straight 
lines. The corresponding deviation for geodesies is 
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which gives a more accurate approximation to geodesies than 
great ellipses. However, the new definition of the spheroidal 
gnomonic projection, Eq. (l93l l. results in an even smaller er- 
ror, Eq. (|96] l, in estimating geodesies. As an illustration, con- 
sider Fig. [13] in which a gnomonic projection of Europe is 
shown. The two circles are geodesic circles of radii 1000 km 
and 2000 km. If the geodesic between two points within one 
of these circles is estimated by using a straight line on this 
figure, the maximum deviation from the true geodesic will be 
about 1.7 m and 28 m, respectively. The maximum changes in 
the end azimuths are 1.1" and 8.6" and the maximum errors 
in the lengths are only 5.4 /im and 730 /im. 

At one time, the gnomonic projection was useful for de- 
termining geodesies graphically. However, the ability to de- 
termine geodesies paths computationally renders such use of 
the projection an anachronism. Nevertheless, the projection 
can be a used within an algorithm to solve some triangulation 
problems. For example, consider a variant of the triangula- 
tion problem 2 of Sect.[TTl determine the point of intersection 
of two geodesies between A and B and between C and D. 



in this case, the method relies on the preservation of azimuths 
about the center point. 

Another application of the gnomonic projection is in solv- 
ing for region intersections, unions, etc. For example the in- 
tersection of two polygons can be determined by projecting 
the polygons to planar polygons with the gnomonic projec- 
tion about some suitable center. Any place were the edges 
of the polygons intersect or nearly intersect in the projected 
space is a candidate for an intersection on the ellipsoid which 
can be found exactly using the techniques given above. The 
inequality (|96] | can be used to define how close to intersec- 
tion the edges must be in projection space be candidates for 
intersection on the ellipsoid. 

The methods described here suffer from the drawback that 
the gnomonic projection can be used to project only about one 
half of the ellipsoid about a given center. This is unlikely 
to be a serious limitation in practice and can, of course, be 
eliminated by partitioning a problem covering a large area into 
a few smaller sub-problems. 



14. MARITIME BOUNDARIES 

Maritime boundaries are defined to be a fixed distance from 
the coast of a state or, in the case of adjacent sta tes or opposite 
states, as the "median line" between the states (iTALOSll2006l 
Chaps. 5-6). In the application of these rule, distances are 
defined as the geodesic distance on a reference ellipsoid to 
the nearest point of a state and the extent of a state is defined 
either by points on the low water mark or straight lines closing 
off bays or joining islands to the mainland. 

For median lines, several cases can then be enumerated 
(iTALOS. ,2006t, Chap. 6): the median is determined by two 
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FIG. 14 The median line (shown as a heavy line) between Britain 
and her North Sea neighbors. Light lines connect the median line 
to the controlling boundary points. The Cassini-Soldner projection 
is used with the central meridian (shown as a dashed line) equal to 
2°20'15"E (the l ongitude of Paris). The data for the coast lines is 
taken from GMT jWessel and SmithlHoiQ) at "intermediate" resolu- 
tion. 



points, a point and a line, or two lines. In the first case, the 
boundary is analogous to the perpendicular bisector in plane 
geometry. It may be constructed by determining the midpoint 
of the geodesic joining the two points and then marking off 
successive points either side of the geodesic by solving the 2- 
distance triangulation problem (problem 1 in Sect. [TTl using 
increasing distances. This continues until some other coastal 
point becomes closer, at which point the median line changes 
direction and continues as the perpendicular bisector of a new 
pair of points. Such turning points are called "tri-points" and 
are equidistant from two points of one state and two point of 
the other; such a point is the center of the geodesic circle cir- 
cumscribing the triangle formed by the three points. 

1 consider first the problem of determining a tri-point O 
given the three coastal points A, B, and C. The solution is 
an iterative one which is conveniently described in terms of 
the azimuthal equidistant projection. Make an initial guess 
O'"' for the position of the tri-point. Map A, B, and C to 
the azimuthal equidistant projection with 0^°^ as the center 
and denote their positions in this projection as a, b, and c. 
Compute the center o of the circle circumscribing the triangle 
formed by these three points. 



Project o to geographic coordinates O*^^^ and use this as the 
new center of projection. Repeat these steps until conver- 
gence. This process converges to the required tri-point be- 
cause of the equidistant property of the projection and it con- 
verges rapidly because of the projection is azimuthal. If the 
points are sufficiently distant (in other words if the center of 
the projection is sufficiently close to the tri-point), Eq. 
can be replaced by 



o = 



((a — h)c + (6 — c)a + (c — a)b) x z 
(a — b) X (b — c) • z 



(98) 



Figure [14] shows the result of using this method to determine 
the median line separating Britain and her North Sea neigh- 
bors. 

With slight modifications this procedure can be applied if 
any of the coast points are replaced by lines. For example, 
assume that A is replaced by a line and let A now denote the 
point on the line closest to O. At the same time as picking 
O'"', provide an estimate A^^^ for the position of A. (These 
can be the result of solving problem 5 of Sect. [TT] however, 
it's more efficient to interleave this solution with the iteration 
for the tii -point.) Transform A^^\ B, and C to the azimuthal 
projection and obtain O'^^ using either Eq. ( |97] l or ( |98l l. Also 
update the estimates for A to A'^^ using one step of Newton's 
method with Eq. (|89] l. In this update, use the computed angle 
between the line and the geodesic from O'"^ to A'"' corrected 
for the anticipated change due to o; this is o x a • z divided by 
the reduced length of the geodesic from 

0(0) to Repeat 

these steps until convergence. 

The median line between two points A and B (or Unes) can 
be found similarly. In this case, introduce a third point G (or 
line), use the distance to C as a control variable, and determine 
the point which is equidistant from A and B and a distance cq 
from C. By adjusting co a set of regularly spaced points along 
the median line can be found. The algorithm is similar to the 
solving for the tii-point using Eq. ( |98l l; however, the updated 
position of the median point in the projected space is 

((c-co)(a-b)-(a-6)c) xz ^^^^ 



(a-b) 



X c • z 



o = 



{a?{h - c) + lP{c - a) + c^{ai - b)) x z 
2(a- b) X (b - c) • z ' 



(97) 



Boundaries which are a fixed distance from a state, typ- 
ically 12 NM for territorial seas or 200 NM for exclusive 
economic zones, can be found using the same machinery 
dTALOSl |2006l, Chap. 5). If the coast is defined by a set of 
points, the boundary is a set of circular arcs which meet at 
points which are equidistant from two coastal points. The 
general problem is to determine the point with is a distance 
ao from A and 6o from B; is a coast point or a point on 
a coastal line and ao = 12 NM; B is another such point 
when determining where two circular arcs meet in which case 
6o = ao. or is a control point in which case 6o measures off 
the distance along a circular arc. This is just problem 1 of 
Sect, [m however, it can also be solving using the azimuthal 
equidistant projection in a similar manner to finding the me- 
dian line. The formula for updating the boundary point in this 
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((a — ao)h — (b — 5o)a) x z 

o = ---^ . (100) 

a X b • z 

Figure[T4]was obtained by exhaustively computing the dis- 
tances to all the coastal points at each point along the median 
line. This is reasonably fast for small data sets, N < 1000 
points. For larger data sets, the points should be organized 
in such a way that the distance to the closest point can be 
computed quickly. For example, the points can be orga- 
nized into quad trees dFin kel and Bentley, 1974) where every 
node in the tree can be characterized by a center and a ra- 
dius r. The geodesic triangle inequality can be used to give 
bounds on the distance s from a point P to any point within 
the node in terms of the distance sq to its center, namely 
max(0, So — r) < s < sq + r. Once the quad tree has been 
constructed, which takes 0{N log N) operations, computing 
the closest distance takes only 0(log N) operations. 



where S12 is the area of the geodesic quadrilateral and inte- 
gral is over the geodesic line (so that cf) is implicitly a function 
of A). Convert this to an integral over the spherical arc length 
a using a similar technique to that used in deriving Eq. (l23T l. 
Sjoberg chooses c — b; however, this leads to a singular in- 
tegrand when the geodesies pass over a pole. (In addition, he 
expresses the integral in terms of the latitude which leads to 
greater errors in its numerical evaluation.) In contrast, 1 define 



2 



2 b^ tanh ^ 1 



(105) 



which is, in effect, the choice that Danielsen makes and which 
leads to a non-singular integrand. The quantity Rg is the au- 
thalic radius, the radius of the sphere with the same area as the 
ellipsoid. Expressing in terms of a gives 

Si2 = Sia2)~SicTi), (106) 
S{a) — R^a + e^a^ cos ao sin 00/4(0-), (107) 



15. AREAS OF A GEODESIC POLYGON 

The last geodesic problem 1 consider is the computation of 
the area of a geo desic polygon. Here, I extend the method of 
lDanielsenl ( ll989h to higher order so that the result is accurate 
to round-off, and I recast his series into a simple trigonometric 
sum which is amenable t o Clenshaw sum mation. In formulat- 
ing the problem, 1 follow [Sioberd (|2006|) . 

The area of an arbitrary region on the ellipsoid is given by 



T 



dT, 



(101) 



where dT = cos (p dcj) dX / K is an element of area and K is 
the Gauss ian curvature. C ompa re this w i th th e Gauss-Bonnet 
theorem (IEisenhartlll940i §34) dBonnetllTsiSL §105) 



KdT, 



(102) 



where F = 27r — 9j is the geodesic excess. This form 
of the theorem applies only for a polygon whose sides are 
geodesies and the sum is over its vertices and 9j is the exterior 
angle at vertex j. Sjoberg combines Eqs. ( llOll l and (1102b to 
give 



— c cos 



(1 — sin^ (j)Y 



d\ 

] cos(j)d<j)d\ (103) 



where c is a constant, and K has been evaluated using Eqs. O 
and ( l25T l. Now apply Eq. ( 11031 ) to the geodesic quadrilateral 
AFHB in Fig.[T]for which V = a2 — a.i and the integration 
over may be performed to give 



Si_2 = (?{oL2 — ai) + 6^ 



A2 



1 



2e sin 



A, v2(l-e2sin^0) 

sin0dA, (104) 



tanh ""^(esin^) (? 



62 



where 



h[cr) = - I ,0 ■ 2 , ;^dcr , (108) 



7r/2 e 



'2 - fc2 sin^ cr' 2 



and 



t{x) = X + \/ X ^ + 1 sinh \px^ 



(109) 



and I have chosen the limits of integration in Eq. ( IIO8I 1 so that 
the mean value of the integral vanishes. Expanding the inte- 
grand in powers of fc2 and e'2 (the same expansion parameters 
as Danielsen uses) and performing the integral gives 



= X] '^4i cos((2; + l)cr) 



(110) 



i=0 
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45 



.10 



792792 ' 



(111) 



I have included terms up to 0{f^) so that the expression for 
S{cr) is accurate to 0{f^). This is consistent with the order 
to which the distance and longitude integrals need to be eval- 
uated to give accuracy to the r ound-off limit for / — 1/150. 
In contrast the s eries given by DanielsenI (1 1989b gives S ac- 
curate to 0{f'^). IClenshawl d 1955b summation can be used to 
perform the (truncated) sum in Eq. (11 lOK with 



ai cos(Z — i)a; = {bi — 62) cos ^x, 

1=1 

where 6; is given by Eq. ( |59] l. 

Summing Eq. ( 1106b for each side of a polygon gives the to- 
tal area of the polygon, provided it does not include a pole. 
If it does, then 27rc'^ should be added to the result. The com- 
bined round-off and truncation error in the evaluation of I4 {a) 
is about 5 x lO^'^ for the WGS84 elUpsoid. However, the 
bigger source of errors is in the computation of the geodesic 
excess, i.e., the first term in Eq. ( 1107) ; the sum of these terms 
gives the area of the spherical polygon and, in Appendix |C] 
1 show how to calculate this accurately. The resulting errors 
in the area of polygon can be estimated using the data for the 
errors in the azimuths given in Sect. |9] Typically, the error is 
approximately a x 15 nm — 0.1 m'^ per vertex; however, it 
may be greater if polygon vertices are very close to a pole or 
if a side is nearly hemispherical. Sometimes the polygon may 
include many thousands of vertices, e.g., determining the area 
of the Japan. In this case an additional source of round-off er- 
ror occurs when summing the separate contributions 5*1 2 frorn 
each edge; this error can be controlled by using Kahanl(ll965b 
summation. 

This method of area computation requires that the sides of 
the polygon be geodesies. If this condition is fulfilled then the 
work of computing the area is proportional to the number of 
sides of the polygon. If the sides are some other sort of lines, 
then additional vertices must be inserted so that the individ- 
ual sides are well approximated by geodesies. Alternatively 
the lines can be mapped to an equal-area projection, such as 
the Albers conic pr ojection (ISnvderill987[ § 14), as suggested 
bv lGilhssenld 1993b . and the area computed in the projected 
space. In either of these cases, the work of computing the area 
will be proportional to the perimeter of the polygon. 



16. CONCLUSIONS 

This paper presents solutions for the direct and inverse geo- 
desic problems which are accurate to close to machine pre- 
cision; the errors are less than 15 mn using double-precision 
arithmetic. The algorithm for the inverse problem always con- 
verges rapidly. The algorithms also give the reduced length 
and geodesic scale; these provide scale factors for geodesic 
projections, allow Newton's method to be used to solve var- 
ious problems in ellipsoidal trigonometry, and enable instru- 
mental errors to be propagated through geodesic calculations. 
I introduce an ellipsoidal generalization of the gnomonic pro- 
jection in which geodesies project to approximately straight 



lines. I discuss the solution of several geodesic problems in- 
volving triangulation and the determination of median l ines. I 
simplify and extend the formulas given by iDanielsenId 19891) 
for the area of geodesic polygons to arbitrary precision. The 
solution of the inverse geodesic problem uses a general so- 
lution for converting from geocentric to geodetic coordinates 
(AppendixlBt. 

Reviewing the list of references, it is remarkable the extent 
to which this paper relies on 19th century authors. Indeed my 
solution of the direct geodesic problems is a straightforward 
extension of that given by Helmert ( 1880) to higher order The 
solution for the inverse problem relies on two relatively mod- 
est innovations, the use of Newton's method to accelerate con- 
vergence and a careful choice of the starting guess to ensure 
convergence; however, the necessary machinery is all avail- 
able in Helmert (1880). My advances, such as they are, rely 
on a few 20th century innovations. The most obvious one 
is the availability of cheap hardware and software for flex- 
ibly carrying out numerical calculations. However, equally 
important for algorithm development are software packages 
for algebraic manipulation and arbitrary precision arithmetic, 
both of which are provided by Maxima (2009) — these facili- 
ties have been available in Maxima since the mid 1970s. 

Computer code implement ing much of th is work is incor- 
porated into GeographicLib (iKarnevi l2010b . This includes 
(a) C++ implementations of the solutions for the direct and 
inverse geodesic problem, (b) methods for computing points 
along a geodesic in terms of distance or spherical arc length, 
(c) computation of the reduced length mi2, the geodesic 
scales M12 and M21, and the area 6*12, (d) implementations 
of the azimuthal equidistant, Cassini-Soldner, and ellipsoidal 
gnomonic projections (all of which return projection scales), 
(e) command-line utilities for solving the main geodesic prob- 
lems, computing geodesic projections, and finding the area of 
a geodesic polygon, (f) Maxima code to generate the series 
Ij (cr) extract the coefficients Aj and Cji, and "write" the C-n- 
code to evaluate the coefficients, (g) the series expansions car- 
ried out to 30th order, and (h) the geodesic test data used in 
Sect. |9] The web page jhttp : //geographiclib.sf . net/geod .html | 
provides quick links to all these resources. In this paper, I have 
tried to document the geodesic algorithms in GeographicLib 
accurately; the source code should be consulted in case of any 
ambiguity. 

I have been questioned on the need for nanometer accuracy 
when geodetic measurements are frequently only accurate to 
about a centimeter I can give four possible answers. (1) Geo- 
desic routines which are accurate to 1 mm, say, can yield sat- 
isfactory results for simple problems. However, more compli- 
cated problems typically require much greater accuracy; for 
example, the two-point equidistant projection may entail the 
solution of ill-conditioned triangles for which millimeter er- 
rors in the geodesic calculation would lead to much larger 
errors in the results. With accurate geodesic routines pack- 
aged as "subroutines", the azimuthal equidistant and Cassini- 
Soldner projections (which are usually expressed as series 
with limited applicability) can be easily and accurately com- 
puted for nearly the whole earth. The need for accuracy 
has has become more pressing with the proliferation of "ge- 
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ographic information systems" which allow users (who may 
be unaware of the pitfalls of error propagation) access to ge- 
ographic data. (2) Even if millimeter errors are tolerable, it 
is frequently important that other properties of geodesies are 
well satisfied, and this is best achieved by ensuring the geo- 
desic calculations are themselves very accurate. An example 
of such a property is the triangle inequality; this implies that 
the shortest path between a point and a geodesic intersects 
the geodesic at right angles and it also ensures the orthogo- 
nality of the polar graticule of the azimuthal equidistant pro- 
jection. (3) Accurate routines may be just as fast as inaccurate 
ones. In particular, the use of Clenshaw summation means that 
there is little penalty to going to 6th order in the expansions 
in Sect. |5] On a 2.66 GHz Intel processor with the g++ com- 
piler, version 4.4.4, solving the direct geodesic problem takes 
0.88 ^s, while the inverse problem takes 2.62 /is (on average). 
Points along a geodesic can be computed at the rate of 0.37 /iS 
(resp. 0.31 /is) per point when the geodesic is parametrized by 
distance (resp. spherical arc length). Thus the time to perform 
a forward and reverse azimuthal equidistant projection (equiv- 
alent to solving the inverse and direct geodesic problems) is 
3.4 fj,s, which is only about 2 times slower than comp uting the 
transv erse Mercator projection using Kriiger's series (Kamevi 
I2OI ih . The object code in geodesic code of GeographicLib 
is substantia l ly long er than an implementation of the method 
of IVincentvl (Il975al) . but, at about 30 kbytes, it is negligible 
compared to the available memory on most computers. (4) Fi- 
nally, it is desirable that a well defined mathematical problem 
have an accurate and complete computational solution; para- 
phrasing Gauss (1903, p. 378) in a letter to Gibers (in 1827, 
on the ellipsoidal corrections to the distribution of the geode- 
sic excess between the angles of a spheroidal triangle), "the 
dignity of science [die Wiirde der Wissenschaft)" requires it. 
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Appendix A: Equations for a geodesic 

Here, I give a derivation of Eqs. ( fT9] l following, f or the 
most part, the presentation of Bessel (1825). Laplacej (1 18291 
Book 1, §8) shows that the path of a geodesic on a surface is 
the same as the motion of a particle constrained to the surface 
but subject to no external forces. For a spheroid of revolution, 
co nservation of an gular momentum gives the relation found 
bv lClairad(ll735h . 

i? sin a = a sin ao , (Al) 

where R is the distance from the axis of revolution (i.e., the 
radius of the circle of latitude), a is the maximum radius of the 
body, a is the azimuth of the geodesic with respect to a merid- 
ian, and ao is the azimuth at the latitude of maximum radius. 
Because R < a, R can be written as a cos /?, where /3 is the 




R=a cosp a 



FIG. 15 The construction for the reduced latitude. The heavy curve 
shows a quarter meridian of the spheroid for which the latitude (f} 
is defined as the angle between the normal and the equator Points 
are transferred from the ellipsoid to the auxiliary sphere, shown as a 
light curve, by preserving the radius of the circle of latitude R. The 
latitude /3 on the auxiliary sphere is the reduced latitude defined by 
R — a cos /?. 

latitude on the auxiliary sphere (see Fig. [15), and Eq. (lAll l 
becomes 

cos /3 sin a = sin ao . ( A2) 

This is the sine rule applied to the angles ao and vr — a in the 
triangle NEP on the auxiliary sphere in Fig.|2]and establishes 
the correspondence with a geodesic on a spheroid with a great 
circle on the auxiliary sphere. It remains to establish the re- 
lations between A and s and their counterparts on the sphere 
u! and a. For a given geodesic on the spheroid, an elementary 
distance ds is related to changes in latitude and longitude by 
(lBesselLll825L Eas. (1)) 

cosQfds = pd0 = — d_R/sin(/i, sinads = -RdA, (A3) 

where p is the meridional radius of curvature. The corre- 
sponding equations on the auxiliary sphere are 

acosadfT = — di?/sin/3, asinada = i?da;. (A4) 

Dividing Eqs. iA3\ by Eqs. ( IA4I ) gives (lBesselLll825L Eqs. (4)) 

1 ds dA sin B 
a da aoj sm (p 

These relations hold for geodesies on any spheroid of revolu- 
tion. Specializing now to an ellipsoid of revolution, paramet- 
rically given hy R — a cos (3 and Z = fosin/3. The slope of 
the meridian ellipse is given by 

dZ b cos /3 cos 4> 
dR a sin /? sin (j) 

This gives the formula for the reduced latitude, Eq. (jS), and 
leads to 

^"^^ f, 2 Tp 

— — - = \/l — e'' cos'^ p = w. 

sin 

Substituting this into Eqs. ( IA5b gives Eqs. ( fT9l ). 
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Appendix B: Transforming geocentric coordinates 



prove the accuracy. The solution proceeds as follows 



IVermeillj (|2002|) presented a closed-form transformation 
from geocentric to geodetic coordinates. However, his so- 
lution does not apply near the center of the earth. Here, I 
remove this restriction and improve the numerical stability of 
the method so that the method is valid everywhere. A key 
equation in Vermeille's method is the same as Eq. (l65T l and this 
method given here can therefore by used t o solve this equa- 
tion. While this paper was being prepared. IVermeillel (1201 ll) 
published an update on his earlier paper which addresses some 
of the same problems. 

As in the main body of this paper, the earth is treated 
as an ellipsoid with equatorial radius a and eccentricity e. 
The geocentric coordinates are represented by {X, Y, Z) and 
the method transforms this to geodetic coordinates (A, 0, h) 
where h is the height measured normally from the surface of 
the ellipsoid. Geocentric coordinates are given in terms of ge- 
ographic coordinates by 

X = {a cos P + h cos 0) cos A, 
Y = {a cos (3 + h cos (f) sin A, 
Z = b sin P -\- h sin 0, 



d = S{S + 2r^), 

u = r + T + /T. 

For d > 0, the sign of the square root in the expression for T 
should match the sign of S* -|- r"^ in order to minimize round- 
off errors; also, the real cube root should be taken. If T = 0, 
then take u — Q. For c? < 0, T is complex, and u is given by 

-0 = ph(-S' - + 1%/^), 

T ^ r exp ^iip, 

u = r(l + 2cosii/')- 

The right-hand side of Eq. (IB 11 1 may now be factored into 2 
quadratic terms in terms of u 



e^K =F ± u), 



(B4) 



where sin/3 — wsm(j), cos/3 = wcos0/Vl — and w is 
given by Eq. ( |20] |. In inverting these equations, the determi- 
nation of A is trivial, 

A = ph{X + iY). 

This reduces the problem to a two-dimensional one, convert- 
ing {R, Z) to ((/), h) where R = \/J0TY^. Vermeille re- 
duces the problem to the solution of an algebraic equation 

At" + 2e^K^ - {x^ + - e^)K^ - 2e^y^K - e^^ = 0, (Bl) 

where 

X — R/a, y = \/ \ — e^Zj a. 



where 



Descartes' rule o f signs shows that for y ^ 0, Eq. (IB II) has 
one positive root (lOlver et al. I l201(i §l.ll(ii)). Similarly for 
X 7^ 0, it has one root satisfying k < — e^. Inside the astroid, 
2;2/3 _|_ y2/3 ^ e'*/'^, there are two additional roots satisfying 
-e^ < n <0. 

The geodetic coordinates are given by substituting the real 
solutions for k into 



(t)^ph{R/{K + e^)+iZ/K), 

h=(i- Vd^ + z^, 



(B2) 
(B3) 



where D = kR/{k + e^). The positive real root gives the 
largest value of h and I call this the "standard solution". 

Equation (IB II) may be solved by standard methods 
(lOlver etal.i \20ld. § l.ll(iii)). Here, I summarize Vermeille's 
solution modifying it to extend its range of validity and to im- 



V ~ \J V? ^ e^y^, 

4 ? 

e y 

V ±u = , for It ^ 0, 

and where the latter equation merely gives a way to avoid 
round-off error in the computation of v ± u. Only the fac- 
tor with the upper signs in Eq. (IB4l i has a positive root given 
by 



V + u 



y^{v + u) + w'^ + ■ 



(B5) 



where 



{v + u)-y'^ 2 

w = e . 

2v 

The number of real roots of Eq. (IBU is determined as follows. 
The condition d ^ is equivalent to x^/"^ + j/^/^ ^ e^^^. 
For d > 0, only the quadratic factor with the upper signs in 
Eq. ( IB4I ) has real roots (satisfying k > and k < — e^, re- 
spectively). For d < 0, both factors have real roots yielding 
the four real roots of Eq. (IB It . 

Equations (IB 21) and ( IB3I ) may become ill-defined if x or y 
vanishes. Solving Eq. dBll ) in the limit a; — > gives 



K^±y, ±e^x/^/e* -y^. (B6) 

Similarly in the limit y ^ 0, Eq. (IB lb yields 

K — —e'^±x, K — zbe^y/V — x^. (B7) 

The only case where this these limiting forms are needed in 
determining the standard solution are for y — and x < e^. 
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Substituting the roots given by the second equation (IB7l l into 
Eqs. (IB2b and ( |B3] ) gives 



= pli(-\/l ^ e^a; ± i\/ — a;^) , 
/i = -6^1 -a;2/e2. 

In the solution given here, I assumed that the ellipsoid is 
oblate. This solution encompasses also the spherical limit, 
e — > 0; the solution becomes +y^. The method may 

also be applied to a prolate ellipsoid, < 0. Substituting 
X = y' , y = x' , K, = k' — in Eq. (IB 11 1 gives 

K —2eK—[x + y — e )K + 2e y k ~ e y =0, 

which transforms the problem for a prolate ellipsoid into an 
equivalent problem for an oblate one. 

In applying the results of this appendix to the inverse geo- 
desic problem, set e = 1 in order to convert Eq. (IBU into 
Eq. I 



Appendix C: Area of a spherical polygon 

The area 5i2, Eq. ( |104| i, includes the term c^{a2 — ai). 
Round-off errors in the evaluation of this term is a potential 
source of error in determining 5*12. In this appendix, I inves- 
tigate ways to compute this term accurately. This term 



a2 — «! 



(CI) 



is of course merely the spherical excess for the quadrilateral 
AFGB in Fig. [T] transferred to the auxiliary sphere. Thus 
Ei2 is the spherical excess for the quadrilateral with vertices 
(/3i,wi), {0,ui), (0,^2), and (^2,i^2)- 

If the geodesic AB is determined by its arc length (T12 and 
its azimuth ai at A then use Eq. ( fT2] i to determine a2 and so 
write 



tan E 12 = 



sin ao cos ao (cos ai — cos (72 ) 
sin'^ ao + cos^ ao cos cti cos a2 



(C2) 



with 



cos (Ti — COS (72 = 

'cosfTi sin(Ti2 



Sm (712 



sm (71 



if cos (712 > 0, 



1 + cos (712 

^cos(7i(l — cos (712) + sin(7i2 sincri, otherwise. 

Here, ao and ai can be determined as described in Sect. |6] 

If the geodesic AB is determined by the latitude and lon- 
gitude of its end points, then, for long arcs, determine ai and 
Q!2 from Eqs. ( l68T l and (|69] l, and substitute these values into 
Eq. dClj ). If, on the other hand, the arc is short, use 



E\2 

tan 

2 



tan i/3i 



1 



1 + tan ^^2 W12 
— — — tan , 

tan ^/?i tan i/32 2 



(C3) 



where, if sin6' and cosO are already known, tan^S may be 
evaluated as sin6'/(l + cos6'),. This relation is the spherical 



generalization of the trapezoidal area; in the limit /?i — > 0, 
(32 — > 0, W12 — > 0, Eq. ( |C3t becomes 

^12 ^ I^12- 

Equation (IC3l l takes on a simpler form if the latitude is ex- 
pressed in terms of the so-called isometric latitude, il> = 

2 tanh^^ tan |/3 = sinh^^ tan /?, namely 

. £^12 ^ , V"! + ■02 , W12 

tan = tanli tan . 

2 2 2 

I obtained Eq. (IC3I) from the formula for the area of a spher- 
ical triangle, E, i n terms of two of t he sides, a and 6, and their 
included angle 7 dTodhuntej. ITSTll §103), 



tan 



E tan ifltan i&sin7 
2 1 + tan ifltan i6cos7' 



by substituting a = ^tt + /3i, & = ^ti" + /32 7 = W12, and 
forming E12 — E — W12. However, it can also be simply 
found by using a formula o f lBessell(fl82l sin. 

q;2 - <3;i sin i(/32 + /3i) ^ a;i2 

tan = r tan , 

2 cosi(/32-/3i) 2 

which, in turn, i s derived from one of Napier's analogies 
(iTodhunteiilTSTTl §52). 

The area for an iV-sided spherical polygon is obtained by 

N 

27m - y^-Ej-i^i, 

i=l 

where n is the number of times the polygon encircles the 
s phere in the ea sterly direction. 

iMillej d 19941) proposed a formula for the area of a spherical 
polygon which made use of L'Huilier's theorem for t he area 
of a spherical triangle in terms of its three sides (Tod hunte^ 
[187 1[ §102). However, any edge of the polygon which is 
nearly aligned with a meridian leads to an ill-conditioned 
triangle which results in about half of the precision of the 
floating-point numbers being lost. 



Appendix D: Geodesies on a prolate ellipsoid 

The focus in the paper has been on oblate ellipsoids. How- 
ever, most of the analysis applies also to prolate ellipsoids 
(/ < 0). In the appendix, I detail those aspects of the problem 
which need to be treated differently in the two cases. 

All the series expansions given in Sect. |5] and the expan- 
sions for S'12 given in Sect. [15] are in terms of /, n, or 
e'2; prolate ellipsoids may be treated easily by allowing these 
quantities to become negative. The method of solving the di- 
rect geodesic problem requires no alteration. The solution of 
inverse problem, on the other hand, is slightly different. From 
Eq. (i23l l, it can be seen that the longitude difference for a geo- 
desic encircling the auxiliary sphere exceeds 2tt. As a con- 
sequence, the shortest geodesic between any two points on 
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region 3 are sufficiently accurate that Newton's method con- 
verges. The geodesic classes in GeographicLib handle prolate 
ellipsoids; however, the algorithms have not been thoroughly 
tested with geodesies on prolate ellipsoids. 

Some of the closed form expressions can be recast into real 
terms for prolate elUpsoids. Thus Eqs. (|36) -(|411 should be 
replaced by 



dn^(u', du' — E{a, k2), 



'0 

1 



F(<7, fca) 



(Dl) 
(D2) 



f Jq 1 — cos^ sn2(u', ^2) 
tan^"'^ (sin ao tan a) 



du' 



1-/ 
/ 



/ sin ao 
G((T, cos^ ao, fc2) 



tan ^(sinaotancr) 
/ sin ao 



(D3) 



where k2 — v— A?, 



am(u2, ^2) 



sn(a;, k) an d dn(.T, k) are Jacobian elliptic functions 

dOlver ef ad I2OIO. §22.2), and G{(j),a'^,k) is defined by 
Eq. (l42b . as before. In this form, the integration constants 
vanish. Finally, Eq. ( 11051 ) becomes 



(D4) 



FIG. 16 (a) Geodesies in antipodal region for a prolate ellipsoid. 
This figure is similar to Fig.|4^, except that / — —1/297. In addition 
the light lines are the continuation of the symmetric set of west-going 
geodesies beyond the meridian A12 — 180°. (b) The dependence of 
A12 on Qi for the near antipodal case with the same value of the 
flattening; compare with Fig.jTJ). 



and Eq. ( 11091 ), which appears in the integral for the geodesic 
area, should be replaced by 



t{-y) = -y + \Jy-^ - Isin ^ 



(D5) 
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